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CHAPTER 1 


INTRODUCTION 

I. Two-Dimensional Thrust Reversing and Vectoring Nozzles 

Two-dimensional nozzles have many advantages for 
fighter aircraft over conventional axisymmetric nozzles. 
Two-dimensional nozzles can be more easily faired to the 
airframe and their use, on twin engine aircraft, eliminates 
the high aft end drag from the "gutter" region between 
axisymmetric nozzles. In addition, two-dimensional nozzles 
can be adapted to thrust reversing and/or thrust vectoring 
with less weight penalty than conventional axisymmetric 
nozzles. 

Thrust reversing and thrust vectoring capabilities are 
likely to be required of future fighter aircraft. The use 
of thrust reversal on landing reduces landing roll and the 
use of thrust vectoring for propulsive lift can reduce both 
takeoff and landing rolls. This short takeoff and landing 
capability is important for operation from bomb damaged 
runways . 

Thrust vectoring and reversing may also be used during 
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flight to increase the aircraft's maneuverability. Thrust 
reversing can modulate the thrust with no spool-up delays, 
and thrust vectoring on aircraft with aft mounted nozzles is 
an effective alternative pitch control. This latter 
capability is important because experimental studies have 
shown considerable loss of elevator effectiveness on some 
configurations when in-flight thrust reversing is used. 

The design of two-dimensional thrust reversing and 
thrust vectoring nozzles requires considerable understanding 
of nozzle internal flow fields. For example, if there were 
large variations in nozzle discharge coefficients during 
thrust reverser deployment it could adversly affect the 
engine performance by causing a mismatch between the choked 
turbine entry area and the effective area of the choked 
main nozzle. Currently, most thrust reversing nozzles are 
designed experimentally. This dissertation 
describes computational procedures developed for the 
analysis of thrust reversing and thrust vectoring nozzles. 

II. Nature of Flow Field 

There are many types of two-dimensional nozzles under 
investigation. The most common of these are convergent 
divergent (CD) nozzles, wedge nozzles, and single expansion 
ramp nozzles (SERN). This investigation deals with thrust 
reversing and thrust vectoring of two-dimensional convergent 
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divergent nozzles (2DCD). 

A typical 2DCD nozzle with thrust reversing and thrust 
vectoring is depicted in Figure 1.1. The gas enters the 
nozzle subsonically and accelerates to transonic speeds as 
it enters one of the three exit ports. The gas then leaves 
the nozzle in a supersonic jet which interacts with the 
ambient flow. Furthermore, there is often a separation 
bubble located on the forward wall of the reverser ports 
near the sharp corners. 

Within two-dimensional nozzles the flow is confined by 
the nozzle sidewalls. There are three-dimensional affects, 
due primarily to the boundary layers on the nozzle 
sidewalls, but these can generally be neglected when the 
nozzle is of high aspect ratio. Once the flow exits the 
nozzle, however, it is no longer confined in the lateral 
direction and the resulting expansion can significantly 
affect the flow within the jet. These three-dimensional 
effects are small near the nozzle exit but will become more 
significant as the distance from the exit increases. 

The principal concern of this investigation is the flow 
within the nozzle. However, for many nozzle geometries the 
flow at the exit plane is subsonic and the nozzle internal 
flow is dependent on the external flow. In these cases the 
external flow field must be calculated. For all cases the 


flow is assumed to be two-dimensional. 


This is a good 



4 


assumption for the internal flow and external flow near the 
nozzle exit. This assumption, however, will cause 
inaccuracies in the external flow far from the nozzle. 
Since the results near the nozzle exit are fairly accurate 
the two-dimensional assumption is believed to be acceptable 
for modelling the effect of external flow on the internal 
flow. 

Since the flow is transonic and contains boundary layer 
separation the effects of both compressibility and 
viscosity must be included in the mathematical model. For 
this reason the two-dimensional compressible Nav ier-Stokes 
equations are solved. These equations are given in Chapter 
2 , Sect ion I . 

III. Related Research 

There have been many experimental studies of 2DCD 
nozzles with thrust reversing and thrust vectoring. Most 
are concerned with gross performance parameters such 
as discharge coefficient and thrust. Two examples are the 
work of Re and Leavitt (Reference 1) and Carson, Capone, and 
Mason (Reference 2). Re and Leavitt studied the static 
internal performance of fully deployed thrust reversing 
nozzles, as well as non-reversing nozzles. Carson, Capone, 
and Mason studied the aero pro pu 1 s i v e characteristics of 
partially and fully deployed thrust reversing nozzles 
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including one partially deployed nozzle with thrust 
vectoring. 

More information concerning the flow field within a 
thrust reversing nozzle is obtained from the measurements of 
Putnam and Strong (Reference 3). They measured static 
pressures along the sidewall, the blocker, and the flap 
(including the forward wall of the reverser port) of a fully 
deployed thrust reversing nozzle. The measurements were 
made with an external ambient Mach number of zero for a 
range of nozzle pressure ratios, from two to eight. 

Several computer programs have been written to solve 
the Na v i er-S tokes equations for flow within conventional 
nozzle configurations. Cline (Reference 4) developed a 
computer program, VNAP, which solved the two-dimensional or 
axisymmetric N a v i e r - S t o k es equations using the 1969 
MacCormack explicit finite difference method. This program 
divides the mesh into two zones so that the mixing of two 
streams can be calculated. Unfortunately, one set of mesh 
lines must always remain vertical, thus limiting the amount 
of geometric turning which may be calculated. VNAP is not 
applicable to thrust reversing nozzles since such nozzles 
invariably have large turning angles. 

Peery and Forester (Reference 16) developed a finite 
volume computer program which also used the 1969 MacCormack 
explicit method. This program utilizes three zones and a 
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generalized nonor thog ona 1 mesh to calculate the mixing of 
three streams. While this code could, in principle, be 
applied to thrust reverser flow fields, no attempt has been 
made to do so. 

Both of the above computer programs solve the Navier- 
Stokes equations using an explicit time marching method. 
Explicit methods are limited by stability to very small time 
steps which are proportional to the smallest mesh spacing 
used. For viscous flow problems the mesh must be refined 
near the wall to resolve the boundary layer profile. As a 
result, explicit methods become very inefficient for viscous 
problems because tens of thousands of time steps are 
required to obtain steady state solutions. Recently, 
implicit methods have been developed which overcome this 
time step limitation. 

Goldberg, Gorski, and Chakravarthy (Reference 5) have 
developed a computer program for axisymmetric afterbody 
flows. This program utilizes a new implicit method which is 
similar to the method presented in this dissertation. The 
program uses an implicit upwind method, with line Gauss 
Seidel relaxation, to solve the Na v ier-Stokes equations on 
a single zone mesh. Since this program does not have 
multiple zone capability it would have limited application 
to thrust reversing nozzles. 

No previous investigators have attempted to solve the 
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N a v ier -S to kes equations for flow with a thrust reversing 
nozzle . 

IV. Contributions of this Research 

A typical thrust reverser flow is very demanding of a 
numerical method. Thrust reverser flows, by nature, involve 
very large turning angles (greater than 90 degrees) and 
rapid expansions. Also, since the thrust reverser is a 
secondary configuration, thrust reversers often contain 
sharp corners where two sections of the baseline forward 
thrust nozzle wall join. Near these sharp corners the rapid 
expansion (with pressure ratios up to five) and large 
turning angles lead to numerical difficulties comparable to 
(or worse than) those encountered while capturing strong 
shock waves. The contributions of this research are 
identifying and overcoming these difficulties. 

For reasons given in the previous section, it was 
decided at the start of this investigation that an implicit 
method was required if the resulting program was to be 
efficient. In theory implicit methods have no time step 
restrictions. However, implicit methods require the 
approximate solution of a very large system of linear 
algebraic equations. Conventional implicit methods 
approximately factor the coefficient matrix for this system 
into the product of two or more simpler matrices whose 
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systems may be efficiently solved using direct methods. 
This approximate factorization leads to an error which is 
small for most flows but significant for thrust reverser 
flows. Our attempts to use various approximate 
factorizations all led to unacceptable time step 
restrictions. The first contribution of this research is 
identifying the approximate factorization as the cause of 
this time step restriction. This problem was overcome by 
solving the system iteratively using line Gauss-Seidel 
relaxat ion . 

The next contribution concerns adapting the solution 
procedure to the complex solution domain found in thrust 
reversing nozzles. Efficient implementation of the line 
Gauss-Seidel relaxation requires computationally simple 
meshes (meshes that may be mapped into a rectangular 
domain). Unfortunately, a single simple mesh cannot easily 
be generated for thrust reverser flow fields like that in 
Figure 1.1. The solution is to break the flow field into as 
many as four zones, each of which is descretized by a simple 
mesh. With a zonal approach such as this, care must be 
taken to insure that the solution procedure is consistent 
across zonal boundaries. This author believes that this is 
the first time a zonal procedure has been utilized with a 
Gauss-Seidel implicit method. 



9 


V. Outline of this Research 

The goal of this research effort was to develop a 
computer program to calculate the viscous compressible flow 
through a two-dimensional con v er gen t- d i v er gen t thrust 
reversing and thrust vectoring nozzle like that shown in 
Figure 1.1. This requires the solution of the compressible 
Nav ier-Stokes equations which are presented in integral 
form in Chapter 2, Section I. These equations are modified 
for time varying meshes in Section II so that the transient 
flows due to the nozzle reconfiguration may be studied. 
Section III presents the algebraic eddy viscosity model used 
to model turbulent flows. 

The solution of the Nav ier - S to kes equations requires 
considerable computer time and an inefficient solution 
procedure may result in such long run times that its use is 
impractical. For this reason an implicit method was chosen 
over explicit methods. Chapter 3 presents the solution 
procedure used. 

Implicit methods generally require that the mesh be 
such that it can be mapped into a rectangular domain. 
Unfortunately, it is very difficult to generate one such 
mesh for the entire flow field in a thrust reversing nozzle 
(Figure 1.1). The alternative is to divide the flow field 
into zones for which simple meshes can be generated. The 
solution procedure is applied to each mesh individually and 
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the zonal solutions are coupled at the zonal interfaces. 
The overall mesh topology and the mesh generation procedure 
for each zone are described in Section II of Chapter 3. 

Section III through V of Chapter 3 present the implicit 
finite volume method used in the interior of the flow field. 
Sections VI and VII discuss the treatment of the boundary 
conditions. Section VIII of Chapter 3 discusses the 
solution of the large system of linear algebraic equations 
resulting from the implicit finite volume method. This is 
one area where this work differs significantly from most 
previous research. Most previous investigations used 
approximate factorization to solve this large system of 
algebraic equations. In Section VIII of Chapter 3, it is 
shown that the standard approximate factorization method 
performs poorly for the thrust reversing nozzle flows. As a 
result, the linear system of equations is solved using a 
line Gauss-Seidel relaxation method which is shown to 
perform dramatically better than the standard approximate 
f act or izat ion . 

In Section IX of Chapter 3, the accuracy and stability 
of the method is analyzed for the model equation. When 
explicit, it is shown to be second-order accurate. When 
fully implicit it is shown to be first-order accurate in 
time, second-order accurate in space, and unconditionally 


stable . 
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In Chapter 4 , the results are presented for five 
nozzle test cases. In Section I the results for a fully 
deployed nozzle are presented and compared with the 
experimental results of Re and Leavitt (Reference 1) and 
Putnam and Strong (Reference 3). The calculated discharge 
coefficient compares very well with the experimental results 
with the largest error being four percent. The calculated 
pressure field also compares well with the experimental 
results, except near the forward wall of the reverser port 
where the pressure is overestimated. This discrepancy 
causes the amount of reverse thrust to be underestimated by 
nearly 15 percent of the ideal thrust when the calculation 
is second order accurate. Numerical dissipation, the 
turbulence model, and three dimensional effects are all 
believed to contribute to this discrepancy. 

In Sections II and III of Chapter 4, results are 
presented for a 50 percent deployed nozzle with both 0° and 
15° of vectoring. Results are compared with the 
experimental data of Carson, Capone and Mason (Reference 2). 
In both cases a coarse mesh was used and the resulting 
numerical dissipation led to an underestimation of the 
discharge coefficient. Both calculations were initially 
first-order accurate which resulted in errors of 21 percent 
in the discharge coefficient. The unvectored case was then 
run second-order accurate and the error was reduced to seven 
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percent. In both cases, the normalized thrust compares 
surprisingly well with experimental data. 

Section IV presents the results for a transient change 
in thrust vectoring angle and Section V the results for a 
transient change in thrust reverser deployment. No 
experimental data is available for these cases, but the 
results look qualitatively correct. 



Upper Reverser Port 
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Figure 1.1 Typical Thrust Reversing and Thrust Vectoring Nozzle 
Flowf ield . 


CHAPTER 2 


Mathematical Model 
I. Navier-Stokes Equation 

The two-dimensional Navier-Stokes equations are a set 
of four coupled nonlinear partial differential equations 
that model the flow of viscous, compressible heat-conducting 
fluids. These equations are derived by application of the 
principles of conservation of mass, momentum, and energy. 
The equations can be written in integral form, 


a_ 

3t 


/■ 


U dV + 


fi 


P»n ds = 


0 , 


where V is an arbitrary volume, S is the surface 
V, and n is the unit vector outward normal to S. 
P = F lx + G Ty 


( 2 . 1 ) 

surrounding 

Also 

( 2 . 2 ) 


where 


P 


U = 


P u 
pv 

_ E _ 


(2.3 a) 
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F = 


pu 


pu + P + T 


XX 


puv + r 


G = 


xy 

(E + p+ r )u+r v + q 
xx xy M x. 

pv 


pvu + r 
r yx 

2 

PV + P + T yy 

r yx u + (e + p + r yy )v + q y 


and 


(2.3 b) 


(2.3 c) 


The variables p, u, v, P, and E are defined below, 
p - density 

u - component of velocity in x-coordinate direction 
v - component of velocity in y-coordinate direction 
P - pressure 

E - volume specific total energy 
The volume specific total energy, E, is related to the mass 

ft 

specific internal energy e, the density p, and the velocity 
by the following equation. 

(2.4) 

The pressure is related to the density and specific internal 
energy by the equations of state, 

P = P ( p ,e) , (2.5) 

which is taken to be the ideal gas law. 

P = P( 7 -1) e (2.6) 


E = p[e + 0 .5 (u 2 + v 2 ) ] 


The gas is also assumed to be calorically perfect so that 
specific internal energy is proportional to the temperature. 
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T, with the constant of proportionality being the specific 
heat at constant volume. 

e = c^T (2.7) 

The specific heat at constant volume, c^, is assumed to be a 
constant equal to 4290.0 ft 2 /s 2 °R, and the ratio of specific 
heats, 7, is assumed to be a constant equal to 1.4. 

To facilitate the development of the solution 
procedure, the contributions to the fluxes from the inviscid 
terms , F^, and the diffusion terms, F^, are considered 
separately. 


F = Fj + 
G = G x + 
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The variables T , T , T , and T are viscous stress 

xx xy yx yy 
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terms which are defined in terms of the derivatives of the 
velocity components. 


T 

XX 

-2 - 

dx 

X V-V 

(2.10 

a) 

'C 

II 

-2M -4— - 

3y 

XV-V 

(2.10 

b) 

II 

T = - 
*yx 
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Ml dy dx 

(2.10 

c) 

v-v = 

3u 3v 

3x + 3y 


(2.10 

d) 


The variables q and q are heat flux terms which are given 

-1 X -3y 3 

in terms of derivatives of temperature. 




(2.11 a) 
(2.11 b) 

The viscosity coef f icients, X and n, are related through 


. 3T 
q = - k -5— 

H y dy 


the bulk viscosity, K . 


3“ ♦ X 


( 2 . 12 ) 

The bulk viscosity is negligible except in the study of the 


structure of shock waves so 


x - - f * 


(2.13) 

The coefficient of viscosity, \i , is related to the 

temperature by Sutherland's empirical formula for air, 

JL = /T ) 3/2 , T, + 198.6 ' <2 ^ 4 ) 

M. 'T. ' T + 198.6' ' U.14; 


where the temperature is expressed in degrees Rankine and 


H o is the reference viscosity at the reference temperature, 
T . The coefficient of thermal conductivity, k, is obtained 
from the expression for the Prandtl number. 


Pr = 


_ Cp/i 


(2.15) 


where Pr is assumed to be a constant equal to 0.72 for 
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laminar flow. 


II. Modified Form for Moving Mesh 

For problems with moving meshes it is convenient to use 
a modified form of these equations. By the Reynolds 
Transport Theorem the first term in equation 2.1 is expanded 
as f ol lows. 

It /" dV - /!f av * utfaoe -" ds 

V VS (2.16) 


Substituting this into the integral form of the Navier- 
Stokes equations (equation 2.1) gives 



* U ''surface 1 •" ds " 0 

s 


(2.17) 


where V £ is the velocity of the surface in the 

surface J 

direction of the unit outer normal, n. Defining a moving 
mesh flux function. 


P 


M 


P + U V 


surface 


(2.18) 


gives 




S 


ds 


= 0 


V 


(2.19) 
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III. Turbulence Model 

For turbulent flow calculations the mass-weighted 
Reynold's averaged N a v i er -S t o k es equations are used 
(Reference 6). These are the standard equations obtained by 
writing the variables in the Na v ier-Stokes as the sum of a 
time averaged quantity and a fluctuating quantity. The 
resulting equations are then averaged and the Boussinesq 
approximation is applied to the resulting Reynolds stresses. 
The result is a set of equations identical in form to the 
equations in section I except that the viscosity and heat 
conduction coefficients become the sum of their laminar 
values with new turbulent coefficients. 

*1 + * t 
k = k^ + k fc 

Here is called the eddy viscosity and k fc is called the 
turbulent conductivity. 

For wall boundary layers, the eddy viscosity is 
obtained from an algebraic two-layer eddy viscosity model 
developed by Baldwin and Lomax (Reference 7). In the 
model, the eddy viscosity, is given by 


( inner 

( n t ) outer 


y' < y' 
y' > y' 


crossover 


crossover 


( 2 . 20 ) 


where y' is the normal distance from the wall and 
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is the smallest value of y* for which ( p ) inner 
crossover t 


= ( /a fc ) outer . 

For the inner layer the Prandtl-Van Driest formulation 
is used. 

( (i t ) inner = p 1^ |o>| (2.21) 

In this equation p is the density, 1 is a length scale, and 
a> is the vorticity. 

( 2 . 22 ) 

The length scale is give by 

1 = k y' [1 - exp(-y + /A + )] (2.23) 

where y + is given in terms of density, shear stress, and 


- g- 

dy dx 


viscosity at the wall. 
+ js/Pw^w y' 


y = 


Pw 


(2.24) 


In the outer region the eddy viscosity is calculated 
with the following formula. 

( outer = K C cp p F wake F KLEB (y') (2.25) 

Here K is the Clauser constant and F T7 , is obtained from 

WAKE 

F WAKE = MIN ^ Y MAX F MAX ' C WK Y MAX U DIFF^ F MAX^* (2.26) 

The terms Y w ,„ and F U11V are obtained from 
MAX MAX 

F(y') = y' |ou| [1 - exp(-y + /A + )] . (2.27) 

where is the maximum value of F(y') in the boundary 

MAX 

layer profile and y MAX is the value of y' at which 

2 

occurs. The term U DIFF 1S simply the square of the maximum 

2 

total velocity, for boundary layer flows. Finally, 

the function F^.-nty') is the Klebanoff i n t erm i 1 1 ency 

L hi i5 
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factor . 


F KLEB (y,) = [1 + 5 * 5 ( — L y B Y )6] 

KLEB y MAX 


-1 


(2.28) 


The constants used in the above formulas are those 
recommended by Baldwin and Lomax (Reference 7). 


A + 

= 

26 

C CP 


1.6 

C KLEB 

= 

0.3 

k 

= 

0.4 

K 

= 

0.0168 


The above turbulence model is implemented for each of 
the interior walls of the nozzle as well as the exterior 
nozzle surfaces. For application to an interior wall the 
testing for and occurs along a column of cells as 
shown in Figure 2.1. This testing occurs from the wall to 
the point half way to the opposite wall (half way in terms 
of number of cells). The spacial derivatives of u and v in 
the expression for the vorticity are calculated using a 
numerical transformation described completely in the 
discussion of the diffusion terms in Chapter 3, Section IV. 

The Baldwin Lomax model is also used for the jets 
emanating from each nozzle exit port. For the jets the 
outer formulation is used with 
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f wake 

and 



(2.29 a) 


f kleb 


(2.29 b) 


Equations 2.29 are utilized with equation 2.25 to obtain the 
eddy viscosities for the jet at each station downstream of 
the nozzle exit. 

In the external flow the results from the jet 
turbulence model must somehow be blended with the results 
from the wall turbulence model. This is done by giving each 
model a domain which it influences as shown in Figure 2.2. 
At the nozzle exit the jet model is used for only those 
cells adjacent to interior cells. For the rest of the cells 
in this column the eddy viscosities are obtained from the 
wall turbulence model. At each successive column of cells 
outward from the nozzle the domain for the jet model widens 
by two cells (one on either side of the jet). After six 
columns the spreading of the jet model increases by four 
cells per column. The wall turbulence model dominates near 
the nozzle (except at the nozzle exits) and the jet model 
dominates far from the nozzle. 

Once the eddy viscosity is known turbulent conductivity 
is obtained from 


"t V pr t 


(2.30) 
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where Pr^_ is the turbulent Prandtl number. As 
Baldwin and Lomax (Reference 7) the turbulent 


recommended by 
Prandtl number 


is taken as 0.9 




I 



Figure 2.2 Domains of Influence for the Jet and Wall Turbulence Models 
in the Exterior Flow. 











CHAPTER 3 


SOLUTION PROCEDURE 
I. Preliminary Comments 

This chapter describes the numerical procedure 
developed for solving the Na v ier-Stokes equations for thrust 
reversing and thrust vectoring nozzle flows. Thrust 
reversing nozzle flows provide difficulties not encountered 
with more conventional nozzles. Typically thrust reversing 
nozzles have sharp corners resulting from the separating of 
the smoothly faired walls present in the cruise 
configuration. Near sharp corners the mesh must be refined 
in two nearly perpendicular directions. This leads to 
unacceptable time step restrictions for conventional 
implicit finite difference methods. The numerical 
procedures presented in this chapter were developed to 
overcome these difficulties. 

Section II describes the mesh topology and mesh 
generation procedure. The mesh is broken into as many as 
four zones to facilitate the application of boundary 
conditions and make efficient use of computer memory. 
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Section III evaluates the time derivative using a 
predictor-corrector method and presents the basic finite 
volume equations. These equations require the evaluation of 
the fluxes through each of the four faces surrounding an 
arbitrary quadrilateral finite volume. Sections IV through 
VI consider the approximate evaluation of these fluxes. 

Section VII presents the application of the boundary 
conditions. Section VIII discusses the solution of the 
large system of linear algebraic equations that result from 
the implicit difference method. It is shown that the 
popular approximate factorization method yields unacceptable 
results for thrust reverser flows and a line relaxation 
procedure is adopted. Finally, Section IX provides an 
accuracy and stability analysis for the solution procedure 
applied to a linear scaler wave equation. 

II. Discretization of the Flow Field 

The flow field under consideration is discretized into 
a large number of finite volume cells. In two dimensions 
these finite volumes are arbitrary quadrilaterals which are 
described completely by the x , y- coo rd i na t es of their four 
corner points. Within the computer program, the mean value 
of the conservative variables, U, within a cell and the 
corner point coordinates are stored in terms of indices i 
and j. Cell i,j, containing U, ., is defined by the x,y- 

* r J 


28 


coordinates with indices i f j, i+l,j, i,j+l, and i+l,j+l. 

On a larger scale, the mesh for the thrust reversing 
and thrust vectoring nozzle program is divided into as many 
as four zones — one zone for each of three possible internal 
flow ports and a zone for the exterior flow as shown in 
Figure 3.1. The mesh is generated such that mesh lines are 
continuous across the zonal interfaces. A single set of i 
and j indices are used to specify cell locations within all 
the internal zones as shown in Figure 3.1 b stacked on top 
of one another. Each must have the same number of columns, 
given by the i index. The mesh cell locations in the 
external flow zone are specified by a separate set of 
indices, ie and je. The outer zone mesh lines must be 
continuous with the internal zone, but the external zone is 
not required to have the same number of columns as the 
internal zones. 

The mesh for each interior zone is generated from cell 
coordinate data input along the lower and upper boundaries 
of the zone and, optionally, also along an arbitrary 
dividing line extending from the inflow boundary to the 
outflow boundary. If the zone has a dividing line, the line 
divides the zone into two regions for mesh generation 
purposes. Otherwise the zone is a single region. An 
example region is shown in Figure 3.2. The mesh is 
generated for this region using a simple algebraic process: 
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1) all corner points for a given i-index within a region 
lie along a line segment connecting the points input for 
that i-index on the upper and lower boundaries, and 

2) the position of an i,j corner point along the line 
segment is given by the relation 



(3.1) 


where 1, L, 1 , and Lj are defined in Figure 3.2. 

The mesh for the external zone is generated somewhat 
differently. The external zone can also be broken into two 
regions as with an internal zone. However, the dividing 
line for the exterior mesh becomes an ie=constant mesh 
line, whereas the dividing line for the interior zone 
becomes a j=constant mesh line. Also, for the interior 
zones, corner point coordinates are given along only three 
of the boundaries of a region where, for the exterior zone, 
corner point coordinates are given for all four boundaries 
of a region. In this case the mesh is generated for each 
region using a modified algebraic process: 

1) all corner points for a given je index within a 
region lie along a line segment connecting the points input 
for that je index along the inner and outer boundaries of 
the region, and 

2) the position of an ie,je corner point along the line 


segment is given by the relation 
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- = (i _ £ — + (- — 

L v A 1 L k a' L 
L u 


(3.2) 


where a is the arclength along the inner boundary of the 
region to point je, A is the total arclength of the inner 
boundary of the region, and 1, L, 1^, L^, l y , and L u are all 
defined in Figure 3.3. 


III. Approximation of the Time Derivative 

The Na v i er-Stokes equations, when applied to a finite 
volume cell, become 


dU . 
1 dt 

i) Vo 1 . 

1 t 

j + D i ( " P * 

S) 

+ Dj (P-S) = 

0 


(3.3) 

where 








D . 
l 

(P- s) = 

P * § i + l/2, j 

- 

P * S i-l/2 , j 



(3.4 a) 

D . 
j 

(P-S) = 

P ’ S i , j+1/2 

- 

P * ^i, j-1/2 



(3.4 b) 

and the 

surface 

vectors are 

always oriented 

i n 

the 

pos i t i ve 

i or j 

directions as shown 

i n 

Figure 3.4. 

If 

t he 

mesh is 


moving the equation is the same. except that P is replaced by 
P M as defined in Chapter 2, Section II. 

The time derivative is approximated with MacCormack's 
predictor-corrector method (Reference 8). 

Pred i ctor 


u n+1 

i/l 


u n 


1 / D 


- (- 


At 

Vo 1 . 


-) [ D. (P- S) P + D. 


(P • S) P ] 


i/3 


(3.5 a) 
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Corrector 
u n+1 = 


0.5 { U n . . + of* 

1 / J 1 9 J 


( voTT~T ) [D i ( P'^) C + D_. (P-s) c n 


1,3 (3.5 b) 

Note the superscripts p and c on the flux terms. These 

indicate that the fluxes may be evaluated differently on the 

predictor and corrector steps. In fact, the way in which 

the fluxes are evaluated in terms of the U. . determines the 

i/3 

type of method that is used. 


IV. Explicit Contribution 

Equations 3.5 do not specify how the fluxes are to be 
evaluated. If the values of the conservative variables were 
known at the surface of each cell (e.g. at i+1/2), the flux 
could be obtained simply from equations (2.2) and (2.3). 
Unfortunately, only the mean value of the solution within 
each cell is known. The flux through a surface must 
therefore be approximated from the mean values of the 
solution within neighboring cells. In general the flux is 
evaluated using both the solution at the current time level, 
nc, and the solution being sought at the new time level, nn. 
This section considers the explicit contribution - the 
contribution from the current time level. 

The inviscid terms and diffusion terms contribute 
substantially different characteristics to the Na v ier-Stokes 
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equations. In the limit as the Reynolds number goes to 
infinity the inviscid terms dominate and the equations are 

hyperbolic in time. In this case information concerning a 

* 

disturbance in the flow field is propagated at a finite rate 
in a manner described by the theory of characteristics. It 
is therefore appropriate to incorporate the theory of 
characteristics into the approximation of the inviscid flux 
through a surface. In the limit as the Reynolds number goes 
to zero the diffusion terms dominate and the equations are 
parabolic in time. The approximation of the surface flux 
due to the diffusion terms should reflect the parabolic 
nature of these terms. This section will consider the 
treatment of the inviscid terms first followed by the 
treatment of the diffusion terms. 

Inviscid Surface Fluxes 

The inviscid flux is evaluated using a second-order 
flux vector splitting which is based, ultimately, on the 
1969 MacCormack method. Three approximations to the surface 
flux are presented here; the 1969 MacCormack method, its 
extension to first-order flux vector splitting, and its 
extension to second-order flux vector splitting. 

1969 MacCormack Method 

In the 1969 MacCormack explicit method (Reference 9) 
the fluxes are evaluated at the latest known time level 
using, alternately, the value of P on either side of the 
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surface. For the i+1/2 and j + 1/2 surfaces (Figure 3.4) 

<f - §) i+ i/2 " f n,j-Vi/2 <3 ' 6 a) 

< p - §) j.i/2- < 3 - 6 b » 

where nc is the current time level, ii is an index which is 
either i or i+1, and jj is an index which is either j or 
j+1. These indices run through a cycle every four steps as 
shown in table 3.1. 

For a subsonic flow, where information can travel in 
either direction, the 1969 MacCormack method violates the 
physical domain of dependence on the predictor and corrector 
steps individually, but satisfies the physical domain of 
dependence collectively. 



1 1 

3D 

np 

nn 

Predictor 

i+1 

j + 1 

n 

n + 1 

Corrector 

i 

j 

n + 1 

n + 1 

Predictor 

i 

j + 1 

n 

n + 1 

Corrector 

i + 1 

j 

n + 1 

n + 1 

Predictor 

i + 1 

j 

n 

n + 1 

Corrector 

i 

j + 1 

n + 1 

n + 1 

Predictor 

i 

j 

n 

n + 1 

Corrector 

i + 1 

j + 1 

n + 1 

n + 1 


Table 3.1 Cycle of Indicies 
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First-Order Flux Vector Splitting 

In general the surface flux should be some nonlinear 
function of the solution on either side of the surface. To 
determine an appropriate approximate function it is 
instructive to look at the eigensystem of the Jacobian 
matrix. First of all, the flux vector is homogeneous of 
degree one in the elements of U. 

(P-S) = [ 1 = AU (3.7) 

The Jacobian can be diagonalized by a similarity 
transformation 

A = T' 1 R -1 S~ 1 A S R T | S | (3.8) 

where S,R, and T are matrices given in appendix A, A is a 
diagonal matrix containing the eigenvalues of the Jacobian 
matrix 

A = 

and u' is the component of velocity in the S direction. 

The eigenvalues of the Jacobian matrix are the speed at 
which information is propogated in the S direction and a 
negative sign indicates propagation of information in the 
direction opposite to S. Figure 3.5 illustrates the 
propagation of information in the S direction for both 


u’ 


u ' 


u ' +c 


u 1 -c 


(3.9) 
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subsonic and supersonic flows. The arrows are the local 
time-dependent characteristics whose slope are given by the 
eigenvalues of A. Following Steger and Warming (Reference 
10) the Jacobian matrix is taken as the sum of two matrices, 
one containing the positive eigenvalues and the other the 
negative eigenvalues. 

A + = T _ 1 R - 1 S -1 A + S R T 
A~ = T - 1 R _ 1 S “ 1 A - S R T 
A + = diag ( A*, A+, A 3 , A*) 

A = drag ( A^ , A^ , A.3 , A^ ) 

N + _ + | ^k | x - - | Ak | 

A k 2 ' A k ' 2 


(3.10 a) 
(3.10 b) 
(3.11 a) 
(3.11 b) 

(3.11 c , d ) 


The flux is then evaluated approximately by multiplying the 
matrix of positive eigenvalues, A + , by the solution to the 
left of the face, and the matrix of negative eigenvalues, 
A , by the solution to the right of the face. 


, - - , . + -,nc - r ,nc 

(P ' S) i+l /2 = A U i,j + A U i+l,j 


(3.12) 


~i_ _ n c 

The jacobians, A and A , are calculated using ^ where 

nc and ii go through the cycle in table 3.1. 


(P- S) 
(P-S) 


FS 

i + 1/2 
FS 

j + 1/2 


if. 


( A 


3 
nc 

i » j j 




3 


(A ) nC . U?* . 

Ur] 1 + lfJ 


u nc . 

1 » J 


(A") ? C 


U 


nc 


jj i/j +1 


(3.13 a) 
(3.13 b) 
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Note that both the positive and negative Jacobians are 
calculated alternately using the values of the solution on 
either side of the surface (i and i+1 for the i+1/2 surface, 
or j and j+1 for the j+1/2 surface). This is in contrast to 
the flux vector splitting of Steger and Warming, (Reference 
10), where * s e valuated using IK j and A i+l/2 1S 
evaluated using U.,, .. As shown in Figure 3.6, the Steger 
and Warming splitting results in five characteristics 
converging on a shock wave and only three characteristics 
converging on a sonic point. Thus the flux at a shock wave 
is overspecified and the flux at a sonic point is 
underspecified. Calculating both Jacobians using the same 
value of U eliminates this problem. 

An additional advantage of this approach is that it 
allows a convenient extension of the 1969 MacCormack scheme 
to flux splitting. The 1969 MacCormack fluxes through the 
i+1/2 face, equations 3.6, may be rewritten as follows. 


.. . . 19 6 9 -nc - r , + . nc . -,nc lr .nc 

n*i: <P-S) i + 1/2 - p i , j * ® i+1/2 = I(A >i,j +(A U i,j 

(3.14 a) 

i i = i + l : (P.8)»« a = * i* 1 , j ‘ § i + 1/2 


i + 1/2 

r( , + ,nc , , -,nc 1 n nc 

= [(A ) i + lfj + (A ) i + lf j]D i + lrj 


(3.14 b) 


Comparing these equations to equations 3.13 reveal that 
the fluxes for first-order flux vector splitting may be 
written as the 1969 MacCormack fluxes plus a second-order 
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smoothing term. 


i i =1 (P-S) 


FS 

i + 1/2 


(P-S) 


1969 
i + 1/2 


+ 


(A )" 


• [U 


n 

i + 1 , j 


ii=i+l (P'S) 


FS 

i + 1/2 


= (P-S) 


1969 
i + 1/2 


- (A + ) 


n+1 
i + 1, j 


[U 


n+1 
i + 1 , j 


u n+1 : 


U i,D ] 

(13.15 


a) 


(13.15 b) 


Clearly, the first-order flux splitting is more 
dissipative than the 1969 MacCormack method. As a result 
the first-order flux splitting is considerably more robust 
than the 1969 MacCormack method, but less accurate. The 
first-order accurate method has other numerical advantages 
when used in an implicit method. These will be discussed 
later . 

Second-Order Flux Vector Splitting 

The accuracy of the flux split method can be improved 
to second-order by using better approximations for U at the 
surface. As shown in Figure 3.7, first-order flux vector 
splitting is equivalent to a zeroth-order extrapolation of 
the solution, U, from the neighboring cell centers to the 
surface. Second-order accuracy is obtained by using a 
linear extrapolation instead (for simplicity the Jacobians 
are evaluated as they were for the first-order flux 
splitting) . 




(D+ »i + l/2 


(3.16 a: 
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/ “ Zi FS2 , ^ ^ \ n c / \ i 

(P ‘ S) j + l/2 = (A ] i, jj (U ] j + 1/2 


(A') nC 

1 i / j j 


(U + ) 


j+1/2 (3.16 


b) 


where 


( ri — ) 

lU i+1/2 

(U + ) nc 
lU ] i+1/2 


.nc 


.nc 


.nc 


UV W . + f , ,, (0V W • - ov~. .) 

1,3 a+ 1/2 i,3 1-1,3 


.nc 


.nc 


“iU.j + <*1/2 ‘"Ul.j - °i + 2,j 


(3.17 a) 
(3.17 b) 


As with the first-order flux vector splitting this is 
obtained as terms that are added onto the 1969 MacCormack 


method . 


.- - . FS2 ,= 5,1969 ,.-,nc r ..nc ..nc , 

li i: (P S ) i + 1/2 ~ (P S) i + 1/2 i , j U i + 1 , j " U i,j 


^ \ o c r ..nc ,i nc n 
* *i.l/2 <A ’ i, j IU i, j - 

.+ ,.-,nc f nc nc , 

^i+1/2 A i,j^i+l,j U i+2,j ] 


(3.18 a) 


ii = i + l: (P S) P “ /2 . (p § )\ll* /2 

/, + .nc rn nc nc -i 

* Vl/2 (A - Vl.j 1 

,+ , nc f „nc ,,nc , 

+ ^i + 1/2 A } i + 1 , j lU i + l, j ~ i + 2 , j 3 


(3.18 b) 


For a mesh with constant spacing 0 i+1/2 and ^i+1/2 
are both 0.5. For a mesh without constant spacing there is 


a choice. The extrapolation can be carried out in physical 

(x,y) space where the extrapolation would depend upon the 

stretching, or it may be carried out in computational (i,j) 

space where the values of <t>~ , , _ and # + , remain 0.5. 

r 1 + 1/2 1 + 1/2 

The latter case is much simpler and is used here. 

Use of the second-order method presented above will 
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result in oscillations near shock waves, rapid expansions, 
and mesh discontinuities. These oscillations lead to 
instabilities near the sharp corner of the thrust reversing 
nozzle when large time steps are taken. A commonly used 
method for eliminating these instabilities is flux limiting. 
With flux limiting the magnitude of the second-order 
contribution is reduced in regions of large gradients. The 
flux limiter used here reduces the coefficients 0 i+1/2 and 
+ by a term proportional to the first and second 

differences of pressure 

*1+1/2 = max ( 0 - 0 ' 0 - 5 ~ CLIM") (3.19 a) 

*1+1/2 = raax ( 0 - 0 ' 0 * 5 “ CLIM + ) (3.19 b) 


where 


CLIM = 


CLIM = 


I p i+1 

- 2 Pi + 

P i-ll 

+ |E>i + 1 

— 

p il 

(p i + l 

+ 2P ^ + 

P i-1 ) 

<p i+ i 


p i } 

l P i.2 

2P i + 1 

+ P il 

. |p w 

- 

p il 

(P i + 2 

+ 2p i.i 

+ Pj) 

(p i+l 

+ 

Pi 5 


(3.20 a) 


(3.20 b) 


This form of flux limiting is equivalent to adding a fourth- 
order smoothing term to the second-order flux split method. 
Diffusion Surface Fluxes 

The viscous stress and heat conduction terms for a cell 
surface, equations 2.9, require the evaluation of spacial 
derivatives of velocity and temperature. These derivatives 
may be obtained from a generalized transformation. 



(3.22) 



Here f and tj are the coordinate directions for a non- 
orthogonal coordinate system with £ running in the i- 
direction and r\ running in the j-direction. The components 
of J^,j can easily be evaluated numerically. For the i + 1/2 
surface , 


tl 

^Ojro 

0 * 25(x i + l, j+ l- X i-l /j+ l +X i + l, 


(3.23 

a) 

8* ... 
8H 

X i+'l,j+l " X i + l,j 


(3.23 

b) 

II 

>l|^ 

0 * 25(y i+l,j+l‘ y i-l,j+l +y i+l, 


(3.23 

c) 

ay _ 

8*7 " 

y i+l , j+1 ” y i+1 , j 


(3.23 

d) 

Also , 





II 

u i+l,j " u i,j 


(3.24 

a) 

8u _ 
an 

0 - 25(u i+l,j+l - u i + l,j-l + U 

i/j+1 ~ U i,j-1 ) 

(3.24 

b) 

The spacial derivatives of v and T 

are evaluated in 

the same 


manner . 


V. Implicit Contribution 


All of the methods presented so far are explicit and 
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are limited to small time steps by the following numerical 
stability criterion: 


At < 


M. + hi + c 

Ax Ay 


+ Ay 2 ] 


172 


( 3 . 25 ) 

When solving the Na v ier-Stokes equations it is necessary to 
have a very fine mesh near walls so that the velocity 
gradients within the boundary layer are adequately resolved. 
Unfortunately, this required that Ax or Ay, and therefore 
At, be very small. In this case the equations are said to be 
stiff and thousands of time steps are required to solve 
interesting problems. 

The difficulty with explicit schemes for Nav ier-Stokes 
solutions arises because explicit schemes require that the 
time step be smaller than the smallest time scale in the 
problem. Look at the mesh near a wall where Ay is the 
direction normal to the wall. Since ^ is much larger than 
^ the equation for the time step becomes 


At < 


Ay 


( 3 . 26 ) 


| v | + c 

Physically this means that the time step must be smaller 
than the time required for an acoustic wave to cross the 
width of the narrowist cell. This time step is unreasonably 
small for two reasons: 

1) The mesh was refined to this degree to resolve 
velocity gradients within the boundary layer, not 
accoustical waves. 
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2) Any inviscid effects this deep in the boundary 
layer should be completely overwhelmed by the viscous 
stresses. The term that limits the stability of the method 
is not even important physically. 

The use of implicit methods can eliminate this time step 
restriction . 

The basic idea of implicit methods is to evaluate the 
fluxes using not only the latest known value of the solution 
as with the explicit methods, but also the unknown solution 
currently being sought. In particular, use a weighted 
average of the fluxes at these two time levels. 

^'^i + i/2 = (1-f i+1/2* i + 1/2 + f i + 1/2 i + 1/2 

(3.27) 

Here fj + i/2 * s t * le ^ e< 3 ree °f implicitness which varies from 
zero (fully explicit) to one (fully implicit) as a function 
of the local CFL number for the surface. 

^ i + 1/2 " MAX 1 0 ' 1 -M > <3 - 28) 

As with the explicit methods the flux terms are still 
evaluated differently between the predictor and corrector 
steps. 

First-Order Flux Splitting 

The fluxes for first-order flux splitting are 
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/r\ 0 \FS > ■. r \ ...nc , -x nc ri nc -» 

p ' s i + l/2 " (1 - f i + l/2> l<A >ii,j u i,j * (A >ii,J u i + l,j> 


+ x nn 


i £ r / * T \ 11 11 r , On . / * — v nil ,, mi 

+ [(A .U. . + (A U., ■ J 

1+1/2 ii,] 1,] li,] 1+1,3 


-,nn 


,nn 


(3.29 a) 


,= ?iFS , > r/ . + .nc ,,nc . , -,nc „nc , 

(p.s ) j+1/2 - (l-f j+1/2 ) [(» + (A irjj^i»j+l 


i -C r , - + v n n ..tin / * \ n n __nn ^ 

+ £ j+l/2 ((A + <A i , jj^i , j+1 


(3.29 b) 


where nc is the current solution, nn is the new solution 
being sought, and ii goes through the cycle in table 
3.1. Substituting these fluxes into equations 3.5 leads to 
a nonlinear algebraic system of equations which must be 
solved on each predictor or corrector step. This is 
impractical, so the system is linearized. For example, take 
the two implicit terms for the i+l/2 surface when ii=i. 


3[ (A + ) . .U . . 

+.nc nc ° l i,D 1,3 
i , . + srr; 


nc 


(A + ) nn .u nn . = (A T )V W .UV~. + — ^ 5U. . + ... 

1,3 1,3 1,3 1,3 3U. . 1,3 

(3.30 a) 


(n ) n nn — (a~ ) 

(A i, j°i+l,j IA ’i,j“i+l,j 


8Ui '5 


nc 


t 3 ltA >i.i u i+l.r 
3D i + l,j 


nc 


«U . 


i + 1 , j 


(3.30 b) 


where 5U • 


i,3 


u nn • - u nc . 

1,3 1,3 


These flux vectors are not homogeneous of degree one so 


but the true Jacobian is expensive to calculate so the 


following is assumed. 


0[ (A ) . . U . .] 

° i,3 i,] 


= (A T ) , 


dU i,3 


if 3 


dt ( a - ) . . u 


ifj i+lfj 


= 0 


30 i,i 


Substituting into equation 3.29 gives 

?\ FS r , + , nc nc , ,.-,nc ri nc , 

(P.S) i+ i /2 . [(A ♦ (A ) ii, j D i+ 1 , i 1 

♦ £ i + i/2 [,s+)n i C i,j 5U i,j * 

(P ’ S1 jtl/2 * [(A+) njj U i?i + IA 

+ £ jti / 2 1 <A+> njj su i.j + 

Substituting th£se fluxes into equations 3.5 leads 

block linear algebraic relationship between 6U . . 

1 , J 

6U in four neighboring cells. 


D . . $U . . . + B . . 5U . . . + A . • 6U . . 

i » 3 i + l f 3 1 , 3 i f 3 + 1 i , 3 1 , 3 

+ C . . 6U . . . + E . . 5U • . . = AU . . 

i,3 1,3-1 1,3 1-1,3 i,3 

where 5U • . is the time difference of U. • obtained 

i,3 i,3 

explicit first order flux split method, and 


D i , j = f i + 1/2 Voli , j 


At_ (A - } nc 

’ ’ n,3 


(3.32 a) 
(3.32 b) 

(3.33 a) 

(3.33 b) 

to a 4x4 
and the 

(3.34) 
from the 

(3.35 a) 
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E . . 

1/3 

B . , 
1/3 

C. . 

i/3 


A . . 

1/3 


- f 

f . 

3 

- f 2 
I + 


At (A + ) nC 

i+1/2 Voli, j ' i i — 1 , j 


At (A~ ) nc 

j + 1/2 Voli , j ; i , j j 
At , +,nc 
j-1/2 Voli , j 1 i / j j-1 
At rf , +>nc 

Voli,j 1 r i + l/2* 7 i i , j 


f i-1/2 (A } 


j -p (, + \ nc 

f j+1/2 A i , j j 


- f . , (A ) nc . . . 

3-1/2 i/33-l 


(3.35 b) 

(3.35 c) 

(3.35 d) 
nc 

i i~l z j 

(3.35 e) 


This relationship is expressed more compactly as 
follows . 

L^ S (5U. .) = AU. . (3.36) 
The subscript I indicates that this is an inviscid 
relationship and the superscript FS indicates that it is 
first-order flux vector splitting. 

When equation 3.36 is applied to each interior cell in 
the finite volume mesh it leads to a large linear system of 
block algebraic equations. To complete this system a block 
algebraic relationship must be provided for each boundary 
cell as well. The boundary cell relationships are developed 
in Section VII of this chapter. The algebraic relationships 
within this system are grouped so that the equation for the 
(i,j) cell is below the equation for the (i,j+l) cell and 
above the equation for the (i,j-l) cell. Thus each column 
of cells is grouped together with equations for the i+1 
column being higher in the grouping than the equation for 
the i column, as shown in Figure 3.8. 

The coefficient matrix for this system of algebraic 
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equations is large, sparse, well structured, and diagonally 
dominant. For example, a 50x50 mesh will result in 
10,000x10,000 coefficient matrices with only 800,000 non- 
zero elements out of 100,000,000 total elements. 
Furthermore, the non-zero elements of the coefficient matrix 
are all contained within 5 bands of 4x4 block matrices, with 
3 of these bands clustered near the main diagonal. The 
principle difficulty with implicit methods is solving this 
large system of equations efficiently. In this 
investigation these systems are solved iteratively using a 
Gauss-Seidel line relaxation method. The motivations for 
using this method, and the details of this method are 
discussed in Section VIII of this chapter. 

Second-Order Flux Spl i tting 

The fluxes for the implicit, second-order flux vector 
splitting method are obtained from equations 3.33 by 
replacing all 50^ by 6uT + 1/2 = «0 lfj * *1.1/2 

“"i.j - < + i / 2 ■ 8 D i + i,j * 

*i+l/2 <8u i+l,3 - 8U i+2,3’ 


i + i/2 ' li, 3 1,3 '1 + 1/2 1,3 i-i/J 

* ,u i + i,j * *m/2 (u i+i,3 - “i + 2,j )))n 

* 8 i+l/2 1 U+) ii, j (8u i ,3 * *1+1/2 <8U i, 3 ‘ 8U i-l 

* (a ‘>U, 3 <8U i + l,j + *i+l/2 (8U i+l,3 - 8U i + 2,i’ 


(3.37 a) 
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(p.s)^/ 2 - «,,, * Vi/2<°i,: - u i, :-i» 

+ (A ) tjjj * *j + l/2 (D l,j + l " °i,jt2 nl 

* f j + l/2 !<»*>“« (SU i , j '*iu/2 lf "i,j - 
+ (a’)"'^ (SU i<j + 1 ♦ *; +1/2 (SD i>j+1 - «0 1>j + 2 ))l 

(3.37 b) 

Substituting these fluxes into equations 3.5 yields a 4x4 
block linear algebraic relationship between 50. and the 5U 
in cells (j+l,j), (i+2,j), (i-l,j), (i-2,j), (i,j+l), 

( i , j + 2 ) , (i/j-1), and (i,j-2). 


D2 . .50. „ . + D. .60. , . + B2 . .80. . - + B. .50. ... 

i»] 1 + 2,] i,] l + l,] i,] i , ] + 2 l,] i,] + l 

+ A. .60. . + C. -50. + C2 • .50. ... + E. -50.,. . 

i,] i,] i,] i,]+l i,] i,]+2 i,] l+l,] 


+ E2 i,: su i+2, j " * u i,: 


(3.38) 


where, AO i . is the time difference of 0 i ^ obtained from the 
explicit second order flux split method, and 


D2 . . 

i ,] 

°i , j 


Voli , j l_f i + l/2 ^i + 1/2 (A } ii,j ] 
Voli, j tf i + l/2 (1 + *i + l/2 )(A } i i , j 


E2 i,r 

E . . 

1 ,] 


+ f i-1/2 *i-l/2 (A } i i-1 , j 3 

V o 1 i , j [ + f i-l/2 *i-l/2 (A ) i-l,j ] 

vfl i , j [_f i-l/2 (1 + *i-l/2 )(A } i i-1 , j 


- f . 


+ * n 


B2 i.j 

B i , j 


Vol i , j C f j + 1/2 ^ j + 1 / 2 ( A } i , j j ] 


i + 1/2 *i + l/2 (A } i i , j 1 


nc 


wrr-i [£ jt i/2« i ♦♦j*i/2"»-)" t : jl 

+ f j-l/2 *j-l/2 (A ] i , j j - 1 ] 


- , nc 


(3.39 a) 

(3.39 b) 
(3.39 c) 

(3.39 d) 
(3.39 e) 


(3.39 f) 
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C2 . 


- At r + f (A + ) nc ] 

Vol i , j +f j-1/2 *j-l/2 lA ; i , j j - 1 J 

At r .c / 1 . f a * ^ n ^ 


: i WT7j '- f 3-l/2 (1 ♦ 

‘ nc 

I . 

rDD 


A . . 

1 f D 


f j+1/2 *j+l/2 (A ) i r jj 1 

1 + VoTT73 [ f i + 1/2 (1 + ^ + l/2 )(A+) ii,j 

- f i -1/2 (1 + ^i-l/2 )(A ' )n i C i-l,j 

+ f i + l/2 (1+ W> (A+) i?jj 

- f , , (1 + 0! , /.) (A ) ? C . . . ] 


(3.39 g) 


(3.39 h) 


(3.39 i) 


This relationship is expressed more compactly as follows. 

l!* S2 (5U. •) = AU • • (3.40) 

1 -*■ f J 1 f J 

The subscript I indicates that this is an inviscid 
relationship and the superscript FS2 indicates that it is 
second-order flux vector splitting. 

The coefficient matrices for this system of algebraic 
equations is large, sparse, and well structured as in the 
first-order case. The example 50x50 mesh yields 
10,000x10,000 coefficient matricies with 1,440,000 nonzero 
elements out of 100,000,000 total elements. The nonzero 
elements of the coefficient matrices are all contained 
within 9 bands of 4x4 block matrices, with 5 of these bands 
clustered near the main diagonal. Unlike the coefficient 
matrices from the first-order flux split method, these 
matrices are not necessarily diagonally dominant. This is 
unfortunate because diagonal dominance is a sufficient 
condition for convergence of the Gauss Seidel line 


relaxation method. 
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Diffusion Terms 

The implicit contribution to the viscous stress and 
heat conduction terms is obtained using a thin-layer 
approximation. To obtain the thin-layer approximation the 
equations are written in Cartesian coordinates aligned with 
the cell face. If^is the coordinate normal to the surface, 
and rj is the cordinate along the surface then the 
contribution to the fluxes from the diffusion terms is 


P-S DIFF 


where 






U'o | + V' 



+ q' 


s 


(3.42) 


u 1 
v' 


=S x u+S y v = 
= -S y u + S x v 


°k = 

nh 


- (A+2M)-|p - X 

/ &u' . 6v' \ 

- - " *' V 

C rp S 

q' - -k 


6T 


6 V 1 
8V 


Velocity normal to surface 
= Velocity tangent to surface 


(3.43 a) 
(3.43 b) 
(3.43 c) 


The thin-layer approximation is obtained by neglecting 
all derivatives in the r? direction. 


a £ = -(X+2 m)-|j- (3.44 a) 

T£v = - (3.44 b) 


The thin-layer diffusion fluxes can then be written. 
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p.s T * L * 


0 


(X+ »>3f S 1 

U &C_ _ . ( , +7u) h^l 


(3.45) 



0 0 0 

0 (X+2p) 0 

0 0 M 

0 u'(X+2p) pv' 


0 

0 

0 

k 


3P 

a£ 

0 u' 

al~ 

av' 

a* 

f 


|S| 

(3.46) 


The first matrix above rotates from the (£, 17 ) coordinate 
system to the (x,y) coordinate system, R ^ . The second is 
the matrix of diffusion coefficients, M. The vector on the 
right contains a set of non- con ser v a t i v e variables in the 
rotated coordinate system, V'. The fluxes may be written in 
a more compact notation. 

P*S T * L * = - R _1 M | S | = - R _1 MR | S | (3.47) 

If£is in the i-direction we can approximate the derivative 
as follows. 


3V . Vy - Vi 2 'Wal (v 

^ (Vol i,5 + Vol i + l,3) i+1 '3 



(3.48) 
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Likewise, if £ is in the j-direction 

, 2|s j^i/2 l (v . . . v ) 

a£ (vol i,j + vol i,j+i ) i,j 


Substituting into equation 3.47 gives 


p.s T ’ L * 

1 b i+l/2 


2|S. 


i + 1/2 


(VOl i + l,j 


; + voi . j (R " lMR) i+i/2 (v i+i, j - v i,j ) 


if] 


(3.50 a) 


p. o T • L • 

P S j+l/2 


( Vol 


2 I ® -i + l/2 I -1 

3 /vol, J (R MR) j + l /2 (V i,j + l- V i,j ) 


i , j+1 


if] 


(3.50 b) 


The purpose of the thin-layer approximation is not to 
get the complete diffusion terms but to provide a simplified 
implicit contribution to the diffusion terms. 


- - ( n T FF ) 

P*S^ + jy 2 = (Explicit Contribution) 


- f . 


2 I S i + l/2 I 


— — : r (R X MR ) . . /0 ( V. . V. .) 

* + Vol, J 1+1/2 i+l,] i f j 


i+l/2(Vol .. 

1 1,3 (3.51 a) 

P *^j + l/2 = Contribution) 

2 I ® -i+ 1/2 I -1 

f j+l/2 (Vol. . + Vol 7 ) (R MR) j+l/2 ( V i, j+l” V i,j } 

1,3 + 1 1,3 (3.51 b) 


The explicit contribution is the full diffusion terms, with 
no thin layer approximation, as given by equation 2.9, 3.23, 

and 3.24. 

The above equations may be written in terms of the 
change in the conservative variables using the 
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(AD *Vl/2 

(AD '>j*l/2 


+ , i 3+1/ + Ot 1/2 - — r (R _ 1 mr) •,!/,«• 

(Vol i,j Vol i , j+1 31/2 3 (3.54 c; 


2f j + 1/2 I S j + 1/2 I ( -1 , 

(V0l i,j + Vol i,j+l ) 3+1/2 j+1 (3.54 d) 


The sum of the flux from second-order flux splitting, 
equations 3.37, and the flux from the diffusion terms, 
equations 3.53, are substituted into equations 3.5. The 
result is a block linear relationship of the same form as 
that discussed in section V. This relationship may be 


written compactly as 
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FS2 


L i + L D <50 i, j ) ■ AU i,j 


(3.55) 


where the explicit diffusion flux is included in AU . .. The 

1 / 3 

diffusion operator is expanded to be a block linear 
relationship between 5U^ ^ and 5U in the four neighboring 
cells . 


L D (5U i, j* - * B i,j 5U i, j+1 + A i, j SO i, j 


(3.56) 


where 


E° . 

B° . 
1/3 

C° . 
i / i 

A° . 
i / 3 


= f . 


Vol i , j (AD * i+1/2 , j 


= f . 


i + 1/2 

f i-1/2 Voli, j (AD } i-1/2, j 
At 


j+1/2 Voli,j ^ AD ^i,j+l/2 


- f . 
3 

At 


j -1/2 Voli, j (AD } i,j-l/2 


Voli , j ff i+1/2 (AD } i+1/2, j f i-1/2 (AD ) 


(3.57 a) 
(3.57 b) 
(3.57 c) 
(3.57 d) 

i-1/2, j 


+ f j + 1/2 (AD ) i, j + 1/2 f j - 1 / 2 (AD ) i,j-l/2 ] 

(3.57 e) 


When the block linear relationship of equation 3.55 is 
applied to each interior cell in the finite volume mesh it 
yields a large linear system of block algebraic relations 
which are ordered in the same manner as the case of second- 
order flux splitting alone. The resulting coefficient 
matrix has nine nonzero block diagonals with five of these 
clustered near the main diagonal. This system of equations 
may be solved in the same manner as the system for the 
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inviscid terms alone. If the thin-layer approximation had 
not been used for the implicit contribution to the diffusion 
terms the resulting system of equations would have been more 
complicated. For the problems considered in this 
investigation, treating only the thin layer terms implicitly 
has proven to yield a stable solution procedure. 


VI. Smoothing Terms 

A second difference of pressure smoother is included to 
allow large time steps during the transient phase of steady 
state calculations. The smoothing is based on simple second 
order smoothing of the conservative variables. 


L (U. .) = VS ... /-> (U • . i • - U. .) - VS. 1/0 (U. .-U. . .) 

s 1,3 1 + 1/2 i+l,3 ir] 1 - 1/2 1,3 i-l,3 

+ VS . (u. . t , - u. . ) - vs. , /0 (u. . - u . ..) 

3 + 1/2 1,3+1 ifD 3 - 1/2 1,3 i,3“l 


(3.58) 


The coefficient for the smoother is proportional to the 
second difference in pressure. 


DDP = 1 P i+3/2 - 2 P i+l/2 * P i-1/2,L (3.59) 

( P i +3/2 + 2 P i +1/2 + P i-l/2 } 

The pressure at a surface is taken to be the average of the 
pressures in the cells adjacent to the surface. This gives 


DDP • 


i +1/2 


I P i + 2 , j P i + 1 , j P i ,1 * P i ~ 1 r j 

<P i + 2 ,i + 3 P i + l,3 * 3 P i,j * P i-l,j' 


(3.60 a) 
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DDP j + l/2 


I P i/j + 2 ~ P i f j + 1 ~ P i,j + 

( P i, j+ 2 + 3P i, j+ l + 3P i fj + P i,j-1> 


VS i + l/2 = max (0r(CS) ( DDP i + 1 / 2 ) " s) 
VS j + l/2 = max ( 0 r (CS) (DDP j + 1/2 ) - S) 


(3.60 b) 


(3.61 a) 
(3.61 b) 


Since the second difference of pressure is proportional 
to the square of the mesh spacing this smoothing is actually 
fourth order. 

The smoother is treated fully implicitly so that the 
U's in equation 3.58 are at the nn time level. Written in 
delta law form the smoothing contribution becomes 


At 


smoothing » vSTT7T [VS i. 1/ 2 ( °W, j ' °l') ) 


,nc 


- VS. 


i-1/2 


+ VS. 


i + 1/2 


- VS 


i-1/2 


+ VS. 


i + 1/2 


- VS . 


i-1/2 


- U T- i,j> 
"^♦i - “Hi' 
,0 “i - 

"“it!,: - «“i,:» 
"“i.: - «“i-i,:» 


+ VS jt 1 /2 
- VS j-l/2 


(s “i,3 - !D i.j-l" 


<5 “i , jtl - *“i, 3 > 

(3.62) 

The first four lines above are the explicit contribution 

from the pressure smoothing and are added into AU. The 

remaining four lines are the implicit contribution which 

form a block algebraic relationship, L ( 5U • .), between the 

S 1 1 3 

five cells with indices ( i , j ) , (i,j+l), (i,j-l), (i+l,j). 



56 


and ( i-1 , j ) . 

When the smoothing in equation 3.62 is added onto 
equation 3.55 the following relationship results. 


L t (5U, .) + L_. ( 5U . •) + L 0 (5U. ■) = AU. 


if 3 


D 


if 3 


3 


i f 3 


(3.63) 


In this relationship the first four lines in equation 3.62 

are included in AU . .. The form of this relation, and the 

1 f 3 

resulting system of algebraic equations, is the same as it 
was before the smoothing terms were added. 


VII. Boundary Conditions 

Explicit 

With the finite volume method the boundaries of a zone 
are placed at cell interfaces and a layer of boundary 
cells surrounds the zone. The only purpose of the boundary 
cells is to satisfy the boundary conditions. For the domain 
considered there are five possible boundary conditions; 
solid wall, plane of symmetry, inflow, outflow, and zonal 
i nter f ace . 

At solid adiabatic walls the pressure and temperature 
gradients are assumed to be zero and a no slip (zero 
velocity at the wall) boundary condition is applied. To 
resolve the boundary layer the mesh must be refined near the 
wall so that the cell nearest the wall is deep within the 
boundary layer. In this case, boundary layer theory shows 
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that pressure gradient normal to the wall is a higher order 
effect and can be neglected. Similarly, for adiabatic walls 
the temperature gradient normal to the wall deep within the 
boundary layer is negligable. These boundary conditions are 
satisfied by setting the pressure and internal energy in 
each boundary cell equal to the pressure and temperature in 
the interior cell adjacent to the boundary. The no-slip 
condition is satisfied by setting the velocity in the 
boundary cell equal and opposite to the velocity in the 
adjacent interior cell. 

The plane of symmetry is a horizontal line at y=0. At 
the plane of symmetry the gradients of pressure, density, 
temperature and x-components of velocity are all zero. 
These conditions are satisfied by setting the pressure, 
density, internal energy, and x-component of velocity in the 
boundary cell equal to their values in the adjacent interior 
cell. The y-component of velocity is zero at the plane of 
symmetry. This condition is satisfied by setting the y- 
component of velocity in the boundary cell to the negative 
of the y-component of velocity in the adjacent interior 
cell . 

The treatment of the inflow boundary conditions is 
guided by the theory of characteristics. A locally one- 
dimensional flow has four characteristic equations with 
slopes u, u + c, u, and u-c. If the flow field is supersonic 
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then all four characteristic equations are propogating 
information in the positive x-direction. In this case all 
data must be specified at the inflow boundary. If the flow 
is subsonic at the inflow boundary, then one of the 
characteristics, the u— c characteristic, has a negative 
slope and it propogates information from the interior 
upstream to the inflow boundary. In this case only three 
items may be specified at the inflow boundary and the fourth 


item must be allowed to vary as the solution progresses. 

For the case of subsonic inflow, the stagnation 
pressure, stagnation temperature, and flow angle are 
specified. These quantities are related to the static 
pressure and static temperature by the following equations: 



1 


7-1 

7+1 



7 

7-1 


( 3 . 64 ) 



7-1 

7+1 



( 3 . 65 ) 


— = Tan(0 xt J = (Constant) 
u IF 


( 3 . 66 ) 


The first two equations above are simply the isentropic 
relations written in terms of the total velocity, V, and the 
speed of sound at a sonic throat, a*. The speed of sound at 
a sonic throat is calculated from the specified stagnation 


I 
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I 

temperature . 

i (a * )2 = 4+1 RT t (3.67) 

Equations 3.64, 3.65, and 3.66 are a system of three 

| equations in four unknowns: p, T, u, and v. To complete the 

system another equation is needed. For all cases considered 
j the flow angle at the inflow boundary is zero so equation 

reduces to v = o and the total velocity, V, is equal to the x- 

i 

I component of velocity, u. 

1 The last equation to close the system is the 

characteristic relation carrying information upstream to the 
r inflow boundary. 

: if -peg - ( U -e»[f -peg] 

} 

I 

This equation is forward differenced. 

[ 

! 5P B - pc 5 u B = (U ~^ )At [Pj ~ P B - pc( Uj - u B )] n (3.69) 

I 

i 

The subscripts I and B indicate the first interior cell and 

I 

! the inflow boundary cell, respectively. The prefix 5 

, indicates the forward in time difference of the variable 

i following it. 

The algebraic equations 3.64, 3.65 and 3.66 can be 
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placed in delta law form by considering incremental changes 
in the variables p, T, and u. 


5p B 


2 Pt u_ 
7 + 1 a* 



(3.70) 


5T B 


7-1 
7+1 a 



(3.71) 


6v = 0 (3.72) 

Equations 3.70, 3.71, and 3.72 are three algebraic 

equations in the three unknowns 5p B , 5 t b , and 8u B . They are 
solved directly for each inflow boundary cell and the 
pressure, temperature, and velocity are updated. 

The delta law formulation of the isentropic 
relations, equations 3.70 and 3.71, simplifies the 
implementation of the inflow boundary. Unfortunately, the 
linearization error allows the stagnation pressure and 
stagnation temperature to vary from the specified values. 
To overcome this problem the static pressure and static 
temperature in each boundary cell is recalculated from 
equations 3.64 and 3.65 using the updated velocity. 

The treatment of outflow boundary conditions is also 
guided by the theory of characteristics. If the flow normal 
to the outflow boundary is supersonic then all 
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characteristics have positive slopes and no information 
propogates upstream from the boundary cell to the interior 
cell. in this case the four locally one dimensional 

characteristic equations are used to update the solution in 
the boundary cell. If the flow normal to the outflow 
boundary is subsonic then the u'-c characteristic propogates 
information upstream from the boundary cell to the interior 
cell. In this case, one item must be specified at the 

boundary cell. 

For subsonic outflow the static pressure is 
specified. The remaining variables in the boundary cell are 
calculated using the three downstream running characteristic 
equations. As with the variables specified at the inflow 
boundary, the requirement of constant pressure at the 
outflow boundary is written in delta law form. 

5P b = 0 (3.73) 

This equation is combined with the three char acter i s ic 
relations . 


8p B + — 5p B = 
c 


u ' At | s | 

' lp B ~ p I + T (P B - P I )] = R 1 


Vol 


I 


(3.74 a) 


5 P b + pc 6 u' b 


( u ' +c)^t s r , , • , , . 

Vol ' ' lP B - P t *PC( 0 ' B - U'j)] = R 2 

(3.74 b) 


u’ At | s 

5 v ' = = 

0 B Vol. 


I 


r v ' - v ' ] = R 

1 B I J k 3 


(3.74 c) 
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The above three equations, along with equation 3.73 are four 
linear algebraic equations in the four unknowns g/t^, fiu' B , 
j v' , and 5p n . This system is solved directly and the 
boundary cell solution is updated. 

At the interzone boundaries the solution in the 
boundary cell is set equal to the solution in the 
corresponding cell in the adjacent zone. 

Implicit 

The explicit boundary conditions described previously 
use the solution at the latest known time level. Implicit 
boundary conditions, however, depend on the solution at the 
unknown time level being sought. This means that the 
solution in the boundary cells is obtained simultaneously 
with the solution in the interior cells. Since the implicit 
boundary conditions are generally linearized in the same way 
as the interior solution procedure the implicit boundary 
conditions will simply be contributions to the large linear 
system of algebraic equations described earlier. The goal 
of this section is to develope implicit versions of the 
boundary conditions described in the previous section and 
write them as block algebraic equations relating the 
solution change in the boundary cell to the change in the 
solution in the adjacent interior cells. 
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The first boundary condition considered is for solid 
walls. For viscous flows the velocity at the wall is zero. 
To impose this boundary condition it is advantageous to 
treat the inviscid fluxes separately from the viscous 
fluxes, in the same manner as these boundary conditions are 
treated explicity. To enforce zero net flux through the 
wall 5U in the boundary cell is set by reflecting the 5U 
from the first interior cell; this in effect is the 
treatment used for free-slip wall boundary conditions. The 
treatment of the viscous fluxes is simply to set the 
velocity components in 6U in the boundary cell to the 
negative of the velocity components in 5U of the first 
interior cell. Separate treatment of the inviscid and 
viscous fluxes at the wall assures that no excessive 
dissipation or spurious numerical fluxes of tangential 
momentum occur. Details of boundary condition treatment are 
described below. 

The inviscid flux treatment for the no-flow through the 
wall boundary condition is satisfied by the equation for the 
boundary cell which relates the change in the solution 
within the boundary cell to the change in the solution in 
the adjacent interior cell. The viscous flux treatment for 
the no-slip boundary condition is implemented by altering 
the viscous Jacobians in the equation for the adjacent 
interior cell. 
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First consider the application of the no-flow boundary 
condition. This condition is satisfied by the following 
relation 


6U ' = EfiU! (3.76) 

O JL 

where 6U' is the change in the conservative variables based 
on velolcities components normal to the boundary, u', and 
tangent to the boundary, v'. The conservative variables 
based on the global coordinate system are obtained from the 
conservative variables based on the local coordinate system 
by multiplying the latter by the inverse of the rotation 
matr ix 


6U = R _1 6U 1 . 


This yields 

5U d = RER _1 5U t . 

D 1 

The E matrix is a reflection matrix: 


E = 


10 0 0 
0-100 
0 0 10 

0 0 0 1 


(3.77) 


(3.78) 


(3.79) 


The coefficient matrices for this block algebraic equation 
are loaded into the global coefficient matrix at the 
position corresponding to the i,j location of the boundary 
cell. For instance, at j=l, i=3 the coefficient matrices 
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are written: 

B 3 1 = - RER _1 (3.80) 

A 3 ^ = I (3.81) 

and all other coefficient matrices for this cell are zero. 

The no-slip boundary condition for the viscous fluxes 
is satisfied by altering the viscous Jacobians. This is 
possible because the dependence of the solution in the 
boundary cell on the solution at the adjacent interior cell 
is known. 


6U b = 

where 


(3.82) 


E NS = 


1 0 

0 -1 

0 0 

0 0 


0 0 

0 0 

-1 0 

0 1 


(3.83) 


Now consider the block algebraic relation for the first 
interior cell and expand it so that the multiplications by 
the viscous and inviscid Jacobians are separate terms. 


. . . + ( A ) 5UJ+ (AD~) 6UJ+ (A + ) 5U B + (AD + ) 5U B + . . . (3.84) 

Use equation 3.82 to write the last term of the above 
equation in terms of 5Uj. 

• . . + ( A~ ) 5U + ( A _ ) 6U j + ( A + ) 5U b + (AD + ) E NS 5 U i + . . . 


(3.85) 
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This approach eliminates the dependence of the implicit 
viscous flux on the change in the solution within the 
boundary cell. The coefficient matrices for the global 
matrix are then modified as follows for j = 2. 


At , A +)n 

, j " j-1/2 Vol . . J i , j j-1 

1 9 J 


C i - " f* 


(3.86 a) 


A i, j 1 


Vof . t *** + f j-l/2 ( (A ) i , j j-1 

1 9 J 

+ (AD')" • . 


_l_ /AH 


M C 


(3.86 b) 


At a plane of symmetry the inviscid terms are treated 
the same as they are at a solid wall (free-slip). The 
viscous terms are treated in a manner similar to that for a 
solid wall, the only difference being the reflection matrix. 
For the plane of symmetry 


5U b = E FS 6U i 

where 


0 0 0 

10 0 
0-10 
0 0 1 


,FS 


(3.87) 


(3.88) 


The plane of symmetry can be located only at j = jl - 1/2 in 
the interior zones or ie = ile -1/2 in the exterior zones. 
Consider only the plane of symmetry at the j = jl - 1/2 


surface . 
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At the interior cell adjacent to the plane of symmetry 
boundary the two modified coefficient matrices for the block 
algebraic relation are 


B 


= f . 


At 


(B )? . 


i » jl-1 j 1-1/2 Vol^ jl _ 1 ' lJ i / j-1 


(3.89 a) 


A i , j 1-1 1 Vol, 


At 


i.ji-i I '" t £ ^- 1/2 




(3.89 b) 


This is completely equivalent to what was done in the no- 
slip case except the reflection matrix is different. 

To treat the inflow boundary condition implicitly the 
differencing of the characteristic relation carrying 
information upstream must be done implicitly. This is done 
by evaluating the spacial differencing using a weighted 
average of the difference at the known time level, nc, and 
the unknown time level, nn. 


5p B - pc6u b = -<l-f 3/2 > IP, - P B - c< u j - V 1 " 0 


B 
- f 


■3/2 
(u-c) At 


Ax 


3/2 Ax 


[ (Pj - p B ) —pc ( u j - u B ) ] 


nn 


(3.90) 


This equation is rewritten. 
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6P b ~ pc6u n = - 


B 


(u-c)At 

Ax 


[ Pj - P B - P c(Uj - u B )] 


nc 


1 « P I - SP B - 0'<* U I - ‘V 1 


(3.91 a) 


I 1 - '‘"b - * ViTT 4 ' '‘"l ■ " c 5 u i ] 

(u-c)At 


Ax 


[p t - p B - ctuj - u B > ) 


nc 


(3.91 b) 


This equation is the implicit equivalent of equation 3.69. 
Provided that Spj and Suj. are known this equation, in 
combination with equations 3.70, 3.71, and 3.72 can be 

solved for 5p B , 5 t 0 , Su B , and Sv B . 

Equation 3.91 b above is used at the inflow boundary 
for the interior zones. For the inflow boundaries in the 
exterior zone a simplified equation is used. This equation 
is obtained by neglecting the second term of equation 3.91b. 


(1 - f 


3/2 


(u-c) At 

Ax 


) I8p b - pc5 u fi ] = 


(u-c) At 

Ax 


[ Pj - p B - pc(Uj 


v ] 


nc 


(3.92) 


The second term of equation 3.91 b is responsible for 
coupling the inflow boundary condition to the interior 
solution during the solution of the linear system. The 
simplified relationship above neglects this coupling. As a 
result, the modified boundary condition can be applied as if 
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it were an explicit boundary condition. Yet it retains the 
unconditional stability of an implicit boundary condition. 
The advantage of this approach is simplicity. 

The outflow boundary conditions are made implicit by 
evaluating R ^ , R^, and R^ in equations 3.74, and using a 
weighted sum of the solution at the known time level nc and 
the unknown time level nn. The result is 


(1 

+ 

\) (6 p B + 

Vp b > - d l U ' , I * V P I> 

c c 

II 

50 

M 3 

(3.93 

a) 

(1 

+ 

d 2 ) (5p B + 

pcfiUg) - d^ ( 8Pj - pc5uj) 

= R n 
R 2 

(3.93 

b) 

(1 

+ 

d l )SV B 

- d l 5v i 

= R n 
R 3 

(3.93 

c) 

(1 

+ 

d 4 ) (sp b - 

pc6u^) - d 4 (fipj - pc5 uj) 

- R n - 
R 4 

supersonic 

/ O Cl 

a \ 


D 1 (3.93 d) 

5p =0 - subsonic (3.93 e) 

B 


where 


At | S | 
d l " Vo — 


(3.94 a) 


d ^ = 


(u J+c) A fc |sj 


Vol 


and 


(uj-c) At |S 
Vol . 


(3.94 b) 


(3.94 c) 


These are four equations relating the change of the 
nonconservative variables in the boundary cell to the change 
of the nonconservative variables within the first interior 
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cell. This bloc-k- linear algehraic relationship is 
incorporated into the global system of equations when the 
external zone is not included. 

When an external zone is included the terms coupling 
the exit boundary solution to the interior solution are 
neglected. The equations for the outflow boundary condition 
then become 


(1 

+ 

diXS^B * Vp b > 

c 

II 


(3.95 

a) 

(1 

+ 

d 2 ) (5p b + pc5u^) 

= R n 
R 2 


(3.95 

b) 

(1 

+ 

<V 5v 'fe 

= R n 
R 3 


(3.95 

c) 

(1 

+ 

d 4 ) (5 P b - pc6u^) 

= R n 
R 4 

- supersonic 

(3.95 

d) 



6 ?b 

= 0 

- subsonic 

(3.95 

e) 


Since these equations are uncoupled from the global system 
they are solved first, as if an explicit boundary condition 
were being used. 

The zonal interface boundary conditions along the lines 
dividing the interior zones are satisfied by simply removing 
from the global coefficient matrix the rows and columns 
corresponding to the zonal interface boundary cells. This 
causes the solution in the cell on one side of the dividing 
line to interact directly with the solution in the cell on 
the other side of the zonal dividing line. 

The zonal boundary conditions for the dividing lines 
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between the interior and exterior zones is handled 
differently. Since the exterior zone uses a different set 
of indices it is more convenient to solve it as a separate 
set of equations. The systems of equations for the interior 
and exterior zones are then approximately solved separately 
and are coupled together through interation by mearly 
setting the U within this boundary cell equal to its latest 
value in the corresponding interior cell of the adjacent 
zone. Since the systems are solved using an interative 
method, this zonal coupling procedure is implemented in a 
manner that is consistent with the overall solution 
procedure. More details concerning interzone boundary 
conditions will be provided in Section VIII. 

VIII. System of Linear Algebraic Equations 

The implicit method presented here requires the solution 
of a large linear system of algebraic equations on both the 
predictor and corrector steps. This system is sparse and 
well structured as shown in Figure 3.81. In particular, 
this system is composed of 9 bands of 4X4 matrices with 5 
of these bands clustered near the main diagonal. The system 
could be solved directly using a Gaussian elimination 
procedure. Unfortunately, this system has a large bandwidth 
and a Gaussian elimination procedure, or variation thereof, 
would fill all of the elements out to the outermost 
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diagonal. For a 50X50 mesh this would require storage for 
8,000,000 floating point numbers, even though only 360,000 
numbers were initially nonzero. Clearly this storage 
requirement is unacceptable. 

The inefficiencies inherent in the direct methods force 
the use of approximate solution procedures. These methods 
generally attempt to approximate the solution of the 
original system by solving a series of simpler systems to 
which direct methods can be applied efficiently. Two 
approximate solution procedures were investigated for the 
linear system of Figure 3.8. In the first, the coefficient 
matrix for the linear system is approximately factored into 
the product of two simpler matices. This approach was 
pioneered by Beam and Warming (Reference 11) and Briley and 
McDonald (Reference 12) for schemes based on central 
differencing. The second procedure approximately solves the 
system using line Gauss-Seidel relaxation. 

The approximate factorization approach is to write the 
coefficient matrix for the system as the product of two 
simpler matrices as shown in Figure 3.9. In this case the 
coefficient matrix in Figure 3.8 is factored into a block 
pentadi agona 1 matrix and a matrix that may be rearranged by 
row and column operations to become a block pe n t ad i a gona 1 
matrix. The B, B2, C, C2, D, D2, E, and E2 elements 
indicated in Figure 3.9 are the same as those defined 
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previously. The AI and AJ matrices are new and are defined 
below. 


AI ifj 1 + 


voTT73 [f in /2 (1 * 


- ^ i -1/2 <1+ *i-l/2 ) <A_) i 1-1 , j 


(3.96 a) 


+ , nc 


i.J - 1 + WITo < f j + l/2 (1 + Vl/2 )(A ’H 

- f j-l/2 11 * *5-l/2> «»'»“!! 


(3.96 b) 


where the A and A Jacobians are defined in equations 3.10 
and the ii and jj indices go through the cycle in table 3.1 
presented before. The two resulting block pentad iagona 1 
systems are solved consecutively using the block 
pentad i agonal solver described in Appendix B. 

There is an error associated with the approximate 
factorization because the product of the two component 
matrices results in a matrix which is different than the 
desired coefficient matrix. The difference of these two 
matrices is the approximate factorization error matrix. 
This matrix is shown in Figure 3.10 for first order flux 
splitting. Each term of the error matrix involves the 
product of an i-direction Jacobian with a j-direction 
Jacobian. Since these Jacobians are proportional to the CFL 
numbers in the i- and j-directions respectively the product 
of these Jacobians is proportional to the product of the CFL 
numbers in the i- and j-directions. In general, then, the 
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approximate factorization method should perform well when 
one or both of the i- or j-direction CFL numbers is small, 
and perform poorly when both the i- and j-direction CFL 
numbers are large. Unfortunately, thrust reversing nozzles 
often contain sharp corners and, near a sharp corner, the 
CFL numbers in both the i- and j-directions are generally 
large. For these cases the approximate factorization method 
has an unacceptable limitation on the time step size. 

The line Gauss-Seidel relaxation method reduces the 
bandwidth of the system by taking either the upper or lower 
off diagonal terms, multiplying them by the solution at the 
previous iteration level, and subtracting the result from 
the right hand side of the equation. This effectively 
applies the contribution of the choosen off diagonal terms 
as a deferred correction. Define IT as the iteration level. 
On the even iteration levels (IT even) the lower off 
diagonal terms are defered and on the odd iteration levels 
(IT odd) the upper off diagonal terms are deferred. If the 
coefficient matrix is such that the line Gauss-Seidel 
relaxation method converges, the approximate solution will 
approach the correct solution as the number of iterations 
becomes large. For first-order flux vector splitting the 
coefficient matrix is diagonally dominant and the line 
Gauss-Seidel relaxation procedure is guaranteed to converge 
(diagonal dominance is a sufficient condition for 
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; convergence of a Gauss-Seidel relaxation procedure). For 

' second-order flux vector splitting the coefficient matrix is 

i not diagonally dominant. Despite this, the line Gauss-Seidel 

relaxation procedure has not diverged for any of the 
problems attempted. 

The implementation of the line Gauss Seidel method is 
best understood by dividing the coefficient matrix into a 
large number of submatrices as shown in Figure 3.8. Each 
row of submatrices is a linear block algebraic equation 
relating the solutions within five adjacent columns of 
cells . 

(D2C)5UC. +2 + (DC)5UC i + 1 + (AC)5UC. + (EOfiUC.^ 

+ (E2C) SUC._ 2 * AUC. (3.97) 

There is one such relationship for each column of cells with 
the submatrix, AC, being a block pentad iagona 1 matrix and 
the submatrices DC, EC, D2C, and E2C being block diagonal 
matrices. Defining IT as the iteration level the line Gauss 
Seidel method is implemented as follows. 

For IT odd: 

1) Set the implicit outflow boundary condition at the 
right boundary (viewed in the computational plane). 

2) Sweep upstream applying the relationship for each 

I T 

column. For the i-column multiply the DC matrix by SUC^ + ^, 
the D2C matrix by 5UcjJ 2 , the EC matrix by SUC^” 1 , and the 
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I T- 1 

E2C matrix by SUC^ 2 . Then subtract the sum of these 

products from the right hand side and solve the block 

I T 

pentadiagonal system for 5UC^ . Repeat this until the left 
boundary is reached. 

For IT even: 

3) Set the implicit inflow boundary condition at the 
left boundary (viewed in the computational plane). 

4) Sweep downstream repeating the process in step 2 

IT-1 IT-1 

except that DC and D2C now multiply 6UC^ + ^ and 8UC^ +2 * and 

IT IT 

EC and E2C now multiply SUc!^ and 5UC^_ 2 . 

5) Repeat the process until the desired iteration 
level has been reached. It was found that only two to four 
iterations were required to achieve stability and that 
overall convergences of the time marching procedure was not 
improved by additional iterations. 

The line Gauss-Seidel relaxation procedure was adopted 
because the time marching scheme based on approximate 
factorization exhibited undesirable instabilities. For 
example, one of the test cases is the symmetric, fully 
deployed thrust reversing nozzle shown in Figure 4.1. The 
key feature of this nozzle is the sharp corner at the 
intersection of the lower flap wall with the forward 
reverser port wall. The initial conditions for the solution 
are stagnation conditions within the nozzle, a low static 
pressure at the exit, and zero velocity everywhere. A 
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physical device that might provide these conditions is a 
diaphram at the nozzle exit plane which is broken at time 
t=0. When the diaphram is broken an expansion wave travels 
upstream through the nozzle until the inflow boundary is 
reached. The inflow boundary conditions then begin pumping 
enough air through the inflow boundary to keep the 
stagnation pressure and temperature at the specified values. 
When approximate factorization is used an instability occurs 
shortly after the expansion wave has passed the sharp 
corner . 

To determine the cause of the instabilities a numerical 
experiment was performed using first-order flux vector 
splitting. The implicit coefficient matrix (for the 
predictor step), based on the transient flow field shown in 
Figures 3.11 through 3.13, was solved using both approximate 
factorization and line Gauss-Seidel relaxation with 2, 3, 4, 
5, 6, and 20 iterations. Figure 3.14 shows the resulting 
solutions as a function of the CFL number for the normalized 
time difference of density, bp/p , in the first cell after 
the corner near the wall. It is clear that the approximate 
factorization procedure gives poor results in this case but 
that the line Gauss-Seidel relaxation gives very good 
results if 3 or more iterations are used. Figure 3.15 
shows similar results for the normalized time difference of 
total energy. While of limited scope, this comparison 
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provides substantial evidence that the Gauss-Seidel 
relaxation procedure is superior for flows with sharp 
corners . 

The time marching scheme based on the line Gauss-Seidel 
relaxation procedure has been successfully applied to this 
fully deployed thrust reversing nozzle. Details of these 
calculations are given in Chapter 4, Section 2, but it is 
significant to note that this method was stable with CF L 
numbers an order of magnitude larger than the CFL number at 
which the approximately factored method was unstable. 

IX. Accuracy and stability Analysis for Model Equation 

The solution procedure presented in the preceeding 

sections is a two-step ex p 1 i c i t- imp 1 ic i t method. This 

method locally varies the degree of implicitness so that the 

procedure is explicit in regions where the explicit 

stability criterion is satisfied, and fully implicit when 

the explicit stability criterion is exceeded significantly. 

In this section the accuracy of this method is analyzed at 

its two extremes (explicit and fully implicit) for the 

following simple model equation (the linear wave equation). 

u. + Xu =0 
t x 

Also in this section, the fully implicit method is shown to 
be unconditionally stable. 
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Accuracy 

First consider the two-step explicit method, with 
spacially second-order upwind differencing, applied to the 
wave equation. 


n+1 n ,3 n 4 n 1 n . 

u . - U, - C(jU. - jUj.j 


“"l 1 - hi * - M‘ n l l - 

where c = 

Ax 


(3.98 a) 
(3.98 b) 


This two-step scheme may be combined to form a single-step 
scheme in the linear case. 


u 


n+1 


n 

i . 
l 


- 

c( l u i 

4 n 

- 2 U i-l 

+ l U i-2> 


c 2 

<K - 

!i u n 

+ ^u" 

- V , 

4 1-3 

2 

4 i-l 

+ 4 U i-2 


♦K-«> 


(3.99) 


To evaluate the accuracy of this method the values of u at 
positions other than (i,n) are approximated using Taylor 
series expansions about point (i,n). Proceeding in this 
manner, it is found that the second term on the right hand 
side of equation 3.99 is a second-order accurate 
approximation to the first partial derivative of u with 
respect to x. 



n 

i 



n 

i-1 


+ 



n 

i-2 


[Ax u - 4 ax 3 u + 0 (Ax 5 ) ] n 
1 x 3 xxx l 


( 3 . 100 ) 


Likewise, the last term is a second-order accurate 
approximation to the second partial derivative of u with 


respect to x. 
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n 

i-1 


+ 



n 

i-2 



n 

i -3 


+ 



n 

i-4 


= [Ax 2 u - ^-Ax 4 u + 0(Ax 5 )]" 

xx 3 xxxx i 


(3.101) 


Equations 3.100 and 3.101 are substituted into equation 

n , 1 

3.99, along with a Taylor series expansion for u ^ . The 

result is 


u t + Xu x 


At , 

- — [u tt 


- X 2 u 


1 + 

xx 3 


. 2 
Ax u 


XXX 


+ H.O.T. 


(3.102) 


The leading term in the truncation error is eliminated by 
differentiating equation 3.102 with respect to t and 
subtracting X times the derivative of equation 3.102 with 
respect to x. 


u„„ - X 2 u = - 
tt xx 2 


[u ttt - Xu ttx - k2 “t„ + 


^“xxx 1 + H -°- T - 
(3.103) 


Substituting into equation 3.102 yields: 


u t + X \ - + - Xu ttx - x2u txx * x3 “xxx : 


+ 7 Ax2 “xxx 


+ H.O.T. 


(3.104) 


The first term of the truncation error in equation 3.104 can 

be eliminated in a similar manner. The reduced equation is 

2 

the wave equation with the leading term being of order Ax . 
u t + \ - I a * 2 “xxx (3 - 105) 
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The two step explicit method with second order upwind 
differencing is therefore second order accurate in both 
space and time. 

Now consider the fully implicit procedure with 
spacially second order upwind differencing. 


n + 1 n r 3 n + 1 4 n+1 , 1 n+1, 

U . 2^., + T U.. 2 ] 


n + 2 n + 1 f 3 n+ 2 

u i - u i - c{ 2 U i 


4 n+2 , 1 n+2, 

2 u i-l + T u i -2 


n+1 l f n , n+2, 

i 2 i i ‘ 


(3.106 a) 


(3.106 b) 
(3.106 c) 


Equation 3.106 may be combined, in the linear case, to 
form a single step scheme. 


n+1 n r 3 n 4n .In. 
1 i - u i * cI 2 u i ' 2 U i -1 + 2“i-2 


r 3 n + 1 

- 2c [ju 
2 

c r - n 


8. 

c_ 

4 


[ 9u 


l 

n + 1 
i 


4 n + 1 
" 2 U i-l 
n 


[ 9 u . - 2 4 u. 


i-1 


^ 1 n+1, 
* 2 u i-2> 

+ 22a "i-2 


n+1 , n+1 

- 24u^_^ + 22u^_2 


8u i-3 + 
o n + l 
- 8u i-3 


n i 
U i-4 
. n + 1 1 
+ u i -4 


(3.107) 


For this equation it is best to expand the Taylor series 
about the point (i,n+l). The results for the four terms 
following the equals sign in equation 3.107 are given below. 


l a i - 7 U i-l + l“i-2 ■ IA, ‘ % - 3itax “tx'f 1 + H -°- T - 


(3.108 a) 


3 n+1 4 n+1 1 n+1 r . „ Ax „ ,n+l 

7 U i - 2 u i-l + 2 u i-2 ■ IAxu x - — "xxx’i 


+ H.O.T. 


(3.108 b) 
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(3.108 c) 
(3.108 d) 

Equations 3.108 are substituted into equation 3.107, along 
with a Taylor series expansion to u”. 

<“t + Xu x> ' % u tt - 6xu tx - x2 “xx» + <3 - 109) 

This equation is reduced in the same way that equation 3.102 
is reduced. The resulting reduced equation is 

u. + Xu = 3AtX 2 u + 0 (At 2 , Ax 2 ) (3.110) 

t X XX 

The fully implicit method with spacially second order upwind 
differencing is second order accurate in space, first order 
accurate in time. 

Stabi 1 i ty 

The stability of the fully implicit scheme is studied 
using a Von Neumann stability analysis. The analysis 
consists of writing the solution as a Fourier series and 
searching for frequencies for which the error grows 
exponentially. 

u ( x , t ) =Zb ( t ) e ikrnx (3.111) 

m 


9 n 

rrU . 

4 1 


24 n , x 22 n 8 t n , In 

4 U i — 1 4 U i-2 " 4 i-3 4 i-4 


= [Ax 2 u xx - 4AtAx 2 


u , ] n+1 + H.O.T. 

txx 1 


9 n + 1 24 n + 1 , 22 n+1 8 n + 1 l, 1 n+ l 

4 U i ' — U i-1 + ~4 U i-2 " 4 U i-3 + 4 U i-4 

r 2 24 , n+1 

= [Ax u xx - T Ax u xxxx ] i 
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Since the model equation is linear, it is sufficient to 
consider an arbitrary term of the Fourier series: 

u m (x,t) = b m (t)e lkmx (3.112) 


This allows the solution at the point (i-l,n) to be 
conveniently written in terms of solution at the. point 
( i ,n) . 


Vi-i't > 


% (x i' tn) e ‘ ls 


(3.113) 


where 0=km x is a frequency parameter. 


Apply the analysis to equations 3.106 


where 


n + 1 


*1 U i 


n + 2 n + 1 n 

1 i = g 2 U i = 9 1 g 2 U i 


n + 1 
u . 

1 


i r - , , n n 

= 2 {1 + g l g 2 } U i = g 3 U i 


g n = g- 


1 + c[ f 


±e~ ifi + ^e -i 0] 
2 e 2 e J 


(3.114 a) 
(3.114 b) 
(3.114 c) 


(3.115) 


For the scheme to be stable the modulus of amplification 
factor, g^, must be less than or equal to one. This is 
equivalent (equations 3.114 c and 3.115) to requiring the 
modulus of g^ to be less than or equal to one. Thus the 
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modulus of the denominator of equation 3.115 must be greater 
than or equal to one. 

||H = | 1 + c [| - ie -i 0 + ie" i2/J ]| > 1 (3.116) 

Expand equation 3.116 in terms of the trigonometric 
functions and multiply by its conjugate. 

2 

| — | = { 1 + c [j - j cos/3 + cos2/3]} 2 

+ c 2 [isin0 - jsin2/3 ] 2 (3.117 a) 

2 

I — I = (Re) 2 + (Im) 2 (3.117 b) 

g l 

2 

The term (Im) is always greater than or equal to zero, 

2 

therefore having (Re) greater than or equal to one is 
sufficient for the stability of the scheme. Expand cos2/3 in 
equation 3.117 a using the double angle formula. 


(Re) = 1 + c[|- - 

4 2 

■jCOS/J + cos /J 


(3.118 

a) 

(Re) = 1 + c( 1 

- cos/3 ) 2 


(3.118 

b) 


From equation 3.118 b it is clear that (Re) is always 
greater than one. The fully implicit method is therefore 
unconditionally stable. 


PHYSICAL PLANE' * COMPUTATIONAL PLANE 


85 


i 

i 



Cm 


Program. 



86 



Total Arclength Along Inner 
Boundary of Region. 

Arlength to Current Line Along 
Inner Boundary of Region. 
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Subsonic Flow (u’>0) 



Supersonic Flow (u*>0) 



Propagation of Information Along Locally One 
Dimensional Time Dependent Characterl sties . 


Figure 3.5 
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Shock Wave - Overspecified 



Sonic Point - Underspecified 



Figure 3.6 Characteristics Resulting from Steger and 
Warming Flux Splitting. 
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First Order Fluxes 



Second Order Fluxes 



Figure 3.7 Extrapolations to the 1+1/2 Surface for First- 
and Second-Order Flux Splitting. 
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A1ISN3Q/ (A1ISN3Q) □ 


figure 3.14 Normalized Time Difference of Density From Approximate 
Solutions of Linear System. 
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A3U3N3/ (A3H3N3) 0 


Figure 3.15 Normalized Time Difference of Energy from Approximate 
Solutions of Linear System. 



CHAPTER 4 


Results and Conclusions 
I. Preliminary Comments 

The computer program described in Chapter 3 has been 
applied to five thrust reversing nozzle configurations. The 
first configuration, described in Section II of this 
chapter, is a fully deployed thrust reversing nozzle with a 
flow turning angle of 130 degrees. The second 
configuration, described in Section III, is a 50% deployed 
thrust reversing nozzle. The third configuration, described 
in Section IV, is a 50% deployed thrust reversing nozzle 
with the forward thrust port vectored downward 15 degrees. 
The final two configurations are transient flow problems. 
In Section V, a calculation of a rapid change in thrust 
vectoring angle is described and in Section VI a calculation 
of a rapid change from partial to full thrust reverser 
deployment is described. 

In Sections II through IV the computed results are 
compared with available experimental data. In all three 
cases the discharge coefficients and normalized thrusts are 
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compared. The discharge coefficient is the mass flow divided 
by the ideal mass flow and the thrust is normalized by the 
ideal thrust. Formulas for the ideal mass flow and ideal 
thrust are given in Appendix C. 

II. Fully Deployed Thrust Reversing Nozzle 

The first test case used during this investigation is 
the fully deployed thrust reverisng nozzle depicted in 
Figure 4.1. This configuration was choosen because 
considerable experimental data is available for it. Re and 
Leavit (Reference 1) have measured the variation of thrust 
and discharge coefficient with nozzle pressure ratio (NPR) 
and Putnam and Strong (Reference 3) have made wall static 
pressure measurements for a series of nozzle pressure 
rat ios . 

The internal flow field for this nozzle was calculated 
using first-order flux splitting for a range of nozzle 
pressure ratios from 2.0 to 7.0. The calculations were 
performed using the 32x27 single zone mesh in Figure 4.2. 
This mesh is divided into two regions for mesh generation 
purposes. The lower region is generated to keep the mesh 
nearly orthogonal to the wall. The mesh is refined along 
the flap and forward reverser port wall to resolve the 
boundary layer. No attempt has been made to resolve the 
boundary layer along the blocker. Since the boundary layer 
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along the blocker is embedded in a highly favorable pressure 
gradient it should be very thin and not significantly effect 
the flow. 

A typical convergence history for these calculations is 
shown in Figure 4.3. This plot shows the variation of mass 
flow error with time step for a nozzle pressure ratio of 
three. The first-order accurate solution required 48 time 
steps to converge. The maximum CFL number for these 
calculations was set to 50,000 for the first four steps, 
100,000 for the next four steps, and 150,000 for the forty 
remaining time steps. 

The flow within this nozzle was also calculated for 
NPR = 3 . 0 using second-order flux vector splitting and the 
mesh shown in Figure 4.4. The internal mesh is the same as 
that in Figure 4.2 and an external mesh was added to 
estimate the effect of the external flow. The second-order 
flux splitting is less robust than first-order flux 
splitting and the CFL number must be set lower. In this 
case the maximum CFL number was set to 2,000. Starting from 
a first-order solution 360 time steps were required to get a 
converged second-order solution. 

The adequacy of the mesh refinement within the 
turbulent boundary layer is presented in Figure 4.5 for the 
second-order calculation with NPR=3.0. This plot gives the 
distribution of the law of the wall coordinate, y + . 


for the 
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mesh points closest to the wall. The variable on the 
abscissa, I, is an index which varies from one at the inflow 
boundary, to 15 just before the sharp corner, to 31 at the 
outflow boundary. Baldwin and Lomax (Reference 7) state 
that a y + of less than two for the mesh point nearest the 
wall, is adequate to resolve the turbulent boundary layer. 
As shown in Figure 4.5, the y + for the mesh point nearest 
the wall is less than two everywhere except at the entrance 
and near the sharp corner. The y + at the entrance is close 
enough to two that the resolution should be sufficient. 
Near the corner the flow is rapidly accelerated and the 
boundary layer becomes very thin. In this case, the 
equilibrium assumption inherent in the Baldwin Lomax model 
is questionable. In any case this region is very small (two 
points with y + significantly greater than two) and it is 
not likely to adversely affect the flow. Overall, the mesh 
refinement near the wall is considered adequate. 

Pressure contours in the region of the exit port for 
nozzle pressure ratios of three, and five are given in 
Figures 4.6 and 4.7, respectively. The pressures, obtained 
using the first-order method, are normalized by their 
respective specified stagnation pressures and all cases had 
the same exit pressure. 

The pressure contours, obtained from the second-order 
method for a nozzle pressure ratio of three, are shown in 
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Figure 4.8. In all three figures the high gradients near 
the corner are evident. The second-order contours. Figure 
4.8, are nearly the same as the first-order contours (Figure 
4.6) for pressures greater than half the stagnation 
pressure. However, there are considerable differences along 
the forward wall of the reverser port with the second-order 
method giving lower pressures than the first-order method. 

Mach number contours in the region of the exit port for 
nozzle pressure ratios of three and five are given in 
Figures 4.9 and 4.10, respectively. Note that the NPR=3.0 
case has a small region of supersonic flow near the exit, 
and the NPR=5.0 case has a much larger region of supersonic 
flow in the reverser port. The Mach number contours, 
obtained from the second-order method with NPR=3.0, are 
shown in Figure 4.11. The second-order method yields higher 
Mach numbers than the first-order method (Figure 4.9). It 
also appears that the boundary layer separation bubble, on 
the forward wall of the reverser port, is larger in the 
second-order case than the first-order case. 

Velocity vectors are also given in the region of the 
exit port. Figure 4.12a shows the velocity vectors for 
NP R = 3 . 0 . A thin layer of separation is evident on the 
forward wall of the reverser port. This separation is shown 
more clearly in Figure 4.12b. Figures 4.13a and 4.13b show 
the velocity vectors for NPR = 5.0. The separation region is 
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smaller than that calculated for NPR=3.0. For comparison, 
the velocity vectors, obtained from the second-order method 
with NPR=3.0, are shown in Figure 4.14a and 4.14b. The 
boundary layer separation obtained from the second-order 
method, Figure 4.14b, is more extensive than that obtained 
from the first-order method. Figure 4.12b. In the first- 
order solution the separation covers 31 percent of the 
forward wall of the reverser port wall and in the second- 
order case it covers 45 percent of the wall. 

The calculated variation of discharge coefficient with 
nozzle pressure ratio is given in Figure 4.15. Also shown 
are the discharge coefficients obtained experimentally by Re 
and Leavit (Reference 1). It is seen that the analysis 
overestimates the discharge coefficient by two to four 
percent. This is to be expected since the effect of the 
sidewall boundary layers is neglected. The side wall 
boundary layers will reduce the effective throat area and 
hence the mass flow. Also shown in Figure 4.15 is the 
second-order calculation at NPR=3.0 (single point). The 
discharge coefficient from the second- order calculation is 
slightly higher than from the first- order method. The 
comparison with the experimental is good for discharge 
coef f icient . 

The calculated variation of thrust with nozzle pressure 
ratio is compared with the experimental results of Re and 
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| Leavitt in Figure 4.16. The first-order results indicate a 

lower amount of reverse thrust than the experimental results 

i 

| with errors ranging from twenty-nine percent of the ideal 

thrust at NPR = 2.0 to eleven percent at NPR = 7.0. Also 
shown in Figure 4.16 is the thrust obtained from the second- 
order method. This thrust compares better with experiment, 
the error being 14.6 percent of the gross thrust compared to 

i 

a 22.0 percent error from the first-order method, 
i The errors in thrust warrant some discussion. The 

error is related to the size of the separation bubble 
calculated on the forward wall of the reverser port. In the 
! calculations the size of the separation bubble seems to be 

underestimated and the pressure on this surface is 

I 

overestimated. This hypothesis is substantiated in part by 
1 the improved results obtained from the second-order method. 

As mentioned earlier, the separation bubble is larger in the 

I 

second-order results than in the first-order results. This 
hypothesis is also supported by the pressure field within 
| the nozzle. 

! Figures 4.17 through 4.19 compare the calculated 

; pressure field with sidewall pressures measured by Putnam 

( 

! and Strong for NPR=3.0. All three figures are for lines, 

parallel to the blocker wall, passing from the plane of 
| symmetry to the exit plane. The line for Figure 4.17 is 

0.203 cm from the blocker wall. At this location the 
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pressures calculated with the second-order method compare 
very well with the measured sidewall pressures. The 
pressures calculated with the first-order method also 
compare well with the experimental results, but not as well 
as the second-order results. Since this line is close to 
the blocker it is an indication that the blocker pressure 
distribution is accurately predicted. Since the 
contribution to the thrust from a surface is approximately 
the integral of the pressure over the surface area, the 
contribution to the thrust from the blocker is also 
accurately predicted. 

The line for Figure 4.18 is 0.838 cm from the blocker 
wall, or approximately half way to the forward wall of the 
reverser port. Again both first-order and second-order 
solutions compare well with the experimental results over 
most of the line. However, the calculated and experimental 
results differ considerably within the exit port. 

The line for Figure 4.19 is 1.437 cm from the blocker 
wall, or nearly to the forward wall of the reverser port. 
Here the sharp corner is identified by the rapid drop in 
pressure. Both the first- and second-order results do well 
up to the corner, but smear out the expansion and miss the 
minimum pressure. The second-order solution is better than 
the first-order solution, missing the minimum pressure by 12 
percent of the stagnation pressure as opposed to a 20 
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percent discrepancy with the* first-order solution. 

The final comparison. Figure 4.20, of calculated 
pressures with experiment is for pressures along the flap, 
rounding the sharp corner, and down the forward wall of the 
reverser port. These pressures are plotted in terms of the 
arclength along the wall from the inflow boundary. The 
experimental results from four lateral stations are also 
plotted. Note that there is considerable variation amongst 
the experimental results, particularly near the exit plane. 
This variation indicates the magnitude of the three 
dimensional effects. It is interesting to note that two 
pressure distributions, both lying within the range of 
experimental pressures presented in Figure 4.20, can yield 
thrusts varying by up to six percent of the ideal thrust. 
The most that can be expected of a two-dimensional solution 
is for it to be within this range. 

Unfortunately, the calculated pressures do not lie 
within the range of pressures measured by Putnam and Strong 
(Reference 3). Both the first- and second-order 
calculations do well before the sharp corner, but predict a 
recompression along the forward wall of the reverser port 
before it is indicated by the experiment results. Thus, 
both the first- and second-order calculations overestimate 
the pressure, and hence force, on the forward wall of the 
reverser port. This accounts for the low calculated reverse 
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thrust compared to the thrust obtained by Re and Leavitt 
(Reference 1). It should be noted that the second-order 
results are considerably better than the first-order 
results. This indicates that the discrepancy between 
calculation and experiment is, in part, due to numerical 
dissipation. The amount of numerical dissipation present in 
a calculation depends on the refinement of the mesh and the 
skewness of the grid, as well as the order and type of 
method used. A coarser or more highly skewed mesh has 
greater numerical dissipation than a finer or nearly 
orthogonal mesh. Unfortunately, the complex geometries 
considered force some degree of skewness to the mesh. The 
refinement, however, is controlled by the user. To obtain 
an accurate solution all features of the flow must be 
resolved. Examination of the velocity vectors in Figure 
4.14b shows that the shear layer at the edge of the 
separation bubble is in a region of coarse mesh. Refining 
the mesh in this region should substantially reduce the 
amount of numerical dissipation in this region, and improve 
the solution. Unfortunately, this calculation was not 
possible due to limited computer resources. 

Another item which might contribute to the error in the 
pressure distribution along the forward wall of the reverser 
port is the forward wall of the reverser port in the 
turbulence model. An overestimation of the eddy viscosity 
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could reduce the size of the separation bubble and cause the 
early recompression seen in Figure 4.20. The Baldwin Lomax 
turbulence model, described in Section III of Chapter 1, has 
been shown to yield satisfactory results for some separated 
flows (Reference 7) but has never been tested on a problem 
as complex as the thrust reversing nozzle. 

Finally, the flow does have significant three- 
dimensional effects as mentioned earlier. The spanwise 
variations in flap pressure distributions (Figure 4.20) are 
probably due to large scale vortices present in the reverser 
port. How these vortices interact with the boundary layer 
separation is not clear. It is possible that the separation 
bubble is enlarged by the vortices. If this is the case an 
accurate prediction of thrust for this geometry can never be 
obtained from a two-dimensional calculation. In any case, 
the two-dimensional calculation does provide good 
quantitative results for the discharge coefficient and 
reasonable qualitative results for the thrust. 

III. Partially Deployed Thrust Reversing Nozzle 

The second test case is the partially deployed thrust 
reversing nozzle shown in Figure 4.21. This nozzle 
geometry, tested by Carson, et.al. (Reference 2), is based 
on a realistic multifunction nozzle design presented by 
Stevens, Thayer, and Fullerton (Reference 13). The 
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experimental results include discharge coefficients and 
thrust for a range of nozzle pressure ratio and ambient 
external Mach number. A calculation is presented for this 
nozzle with a nozzle pressure ratio of two and an ambient 
external Mach number of zero. 

The calculation was performed using the mesh shown in 
Figure 4.21. Due to limited computer resources, this mesh 
is very coarse and none of the wall boundary layers are 
adaquately resolved. 

The resulting discharge coefficient is .89 compared 
to .993 from the experiment of Carson, et.al. There is 
evidence that the mesh is inadaquate for accurately 
calculating the mass flow. The problem was run to 
convergence using the first-order accurate results and then 
switched to second-order accuracy. The first-order result 
for discharge coefficient is only .78 compared to .89 from 
second-order accuracy. If the mesh were sufficiently fine 
the difference would be less dramatic. 

It should be noted that the areas of the reverser ports 

for the model tested by Carson, et.al. are somewhat larger 

than for the nozzle used in the calculation. the combined 

areas of the three exit ports in the experiment were 36.11 
2 2 

cm as opposed to 34.89cm for the calculation (based on two 
nozzles, each with a width of 7.77 cm). Carson, et.al. 
based their ideal mass flow calculation on the throat area 
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2 

of the baseline forward thrust nozzle, 34.24 cm . 
Correcting the calculated discharge coefficient (by simply 
multiplying by the ratio of areas) yields a discharge 
coefficient of 0.92. 

The calculated thrust is found to be slightly higher 
than the measured thrust. The calculated thrust, normalized 
by the ideal thrust based on the measured mass flow, is 
0.253, whereas the measured thrust, normalized by the ideal 
thrust based on the measured mass flow, is 0.23. This 
comparison is suprisingly good considering the coarseness of 
the mesh. 

The pressure contours, Mach number contours, and 
velocity vectors for this nozzle are given in Figures 4.22 
through 4.24. It is clear from Figure 4.23 that the flow 
remains subsonic until it leaves the exit port. Thus, the 
internal flow is influenced by the external flow and the 
external zone is needed. 

IV. Partially Deployed Thrust Reversing Nozzle with 
Vectoring 

The third test case, shown in Figure 4.25, is the same 
50% deployed thrust reversing nozzle discussd in Section III 
with the forward thrust port vectored 15° downward. This 
nozzle was also tested by Carson, et.al., and the discharge 
coefficient and thrust-were obtained for a range of nozzle 
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pressure ratios and two free stream Mach numbers. The 
calculation was made at a nozzle pressure ratio of two with 
an ambient free stream Mach number of zero. 

As with the unvectored case the calculation was 
performed using a coarse mesh. The mesh is composed of four 
zones: three internal zones and an external zone. The 
dimensions of the mesh for all three internal zones is 
21x10. For the external zone the mesh dimensions are 
17x107. 

The calculation was run to convergence first order 
accurate (64 time steps). The resulting discharge 
coefficient is 0.77, which is consistent with the results 
for the unvectored case. The calculated thrust, normalized 
by the ideal thrust, is 0.146 and the normal force, 
normalized by the ideal thrust is 0.065 Note that the ideal 
thrust is based on the measured mass flow (Appendix C). 

The thrust and normal force normalized by the ideal 
thrust based on the calculated mass flow are 0.19 and 0.085. 

The experimental results for this case had considerable 
asymmetry. The thrust was measured for vectoring of 15° 
both upward and downward. With nozzle flow vectored 
downward the normalized thrust was 0.18 whereas it was 0.34 
when vectored upward. The computed normalized thrust is at 
the lower end of this range. Likewise, the computed jet 
turning angle of 24° lies between the experimental results 
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of -5° and 45°. 

The calculated pressure contours, Mach number contours, 
and velocity vectors for this case are shown in Figures 4.26 
through 4.28. It appears, from the Mach contours in Figure 
4.27 that there is no supersonic flow. The results for the 
unvectored nozzle, however, indicate small regions of 
supersonic flow just outside each exit port when second- 
order accuracy is used. It appears that the first-order 
accuracy, in combination with the coarse mesh, has resulted 
in excessive total pressure loss. This seems to be the 
reason for the low discharge coefficient. If this 
calculation was done using second-order accuracy, or the 
mesh was refined, the calculated discharge coefficient would 
be much closer to the measured discharge coefficient. 

V. Transient Change in Thrust Vectoring Angle 

A calculation was made of a transient change in thrust 
vectoring angle. The nozzle was initially unvectored and 
the vectoring angle was changed from 0° to 30° over a period 
of .3 msec. The actual wall velocity followed a cosine 
function of time with the maximum velocity being nearly half 
the speed of sound. This very fast change was necessary for 
transient effects to be observed. The initial solution is 
shown in Figures4.29 through 4.31. The mesh uses three 
internal zones with the lower and upper zones having 
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dimensions 23x16 and the middle zone having dimensions 23x8. 
The nozzle pressure ratio is 5.0. 

The effects of the moving wall are seen by comparing 
the pressure contours during the final stages of the 
reconfiguration. Figure 4.34, with the pressure contours 
after the reconfiguration is complete. Figure 4.37. The 0.2 
contour in Figure 4.34 extends from the exit of the rear 
port half the distance to the corner. In Figure 4.37 this 
contour has relaxed so that it connects to the wall near the 
exit port and intersects the exit plane at a higher 
position. No experimental data is available for this case. 

VI. Transient Change in Thrust Reverser Deployment 

A calculation was made of a transient change in degree 
of thrust reverser deployment. The nozzle was initially 70 
percent deployed with the rear port width being 0.05 ft. 
The rear port was then closed, as the reverser ports were 
opened, over a period of .3 msec. 

The results, shown in Figures 4.38 through 4.46, 
indicate much more dramatic transient affects than the 
transient thrust vectoring case. As shown in the pressure 
contours, (Figure 4.43) the closing rear port compresses the 
air near the port to much greater than the stagnation 
pressure at the inflow boundary. As the port closes the 
pressure becomes large enough that the flow near the 
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entrance to the port reverses direction and the air flows 
back into the nozzle as shown in Figures 4.42a and 4.42b. 
There is no experimental data available for this case. 

VII. Conclusions 

An implicit finite volume program has been developed 
for the calculation of two-dimensional thrust reversing 
nozzle flows. Thrust reversing nozzles typically have sharp 
corners, and the rapid expansion and large turning angles 
near these corners lead to unacceptable time step 
limitations when conventional approximate factorization 
methods are used. In this investigation these limitations 
are overcome by replacing the approximate factorization with 
a line Gauss-Seidel relaxation. This method is implemented 
with a zonal mesh so that flows through complex nozzle 
geometries can be efficiently calculated. 

Results are presented from calculations using both 
first- and second-order fully upwind differencing. In most 
cases, the second-order method compared better with the 
experimental results than the first-order method. The 
second-order method was limitted, by stability, to lower CFL 
numbers than the first-order method. The lower CFL numbers 
resulted in longer run times for the second-order method. 
(Typically, 360 time steps for the second-order method 
compared to 48 time steps for the first-order method). 
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Despite this limitation, the second-order method still 
requires two to three orders of magnitude fewer time steps 
than an explicit method when turbulent boundary layers are 
adequately resolved. 

The comparisons of the calculated results with 
experiment were mixed. For the fully deployed nozzle 
(Section II of this chapter) the calculated discharge 
coefficient compared well with experiment but the amount of 
reverse thrust was underestimated. Conversely, for the two 
partially deployed nozzles (Sections III and IV) the 
calculated thrust compares well with experiment but the 
discharge coefficient is underestimated. The 
underestimation of the discharge coefficient for the 
partially deployed nozzles is undoubtably due to the coarse 
mesh used for the calculation. The underestimation of the 
reverse thrust for the fully deployed nozzle is thought to 
be caused by a combination of insufficient mesh density, 
limitations of the turbulence model, and three-dimensional 
effects in the experiment which cannot be predicted with a 
two-dimensional model (see Section II of this chapter). 

The computer program is robust, efficient, and capable 
of calculating the complex flows within thrust reversing 


nozzles . 




Figure 4.1 Details of 2-D Convergent-Divergent Thrust-Reversing 
Nozzle. Internal width of nozzle is 10.160 cm. All 
Dimensions are in centimeters unless otherwise noted 
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Figure 4.3 


Hass Flow Convergence History for the 
Fully Deployed Thrust Reversing Nozzle 







CORNER 



Evaluation of Grid Refinement Along Flap. Y + Distribution 
for Cells Nearest Wall. 
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Figure 4.6 Pressure Contours for Fully Deployed Thrust Reversing 
Nozzle - NPR = 3.0 - 32x27 Mesh - First-Order. 




Figure 4.7 Pressure Contours for Fully Deployed Thrust Reversing 
Nozzle - NPR = 5.0 32x27 Mesh - First-Order. 
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Figure 4.8 Pressure Contours for Fully Deployed Thrust Reversing 
Nozzle - NPR = 3.0 - 32x27 Mesh - Second-Order. 



umber Contours for Fully Deployed Thrust Reversing 
- NPR = 3.0 - 32x27 Mesh - First-Order. 
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Figure 4.10 Mach Number Contours for Fully Deployed Thrust Reversing 
Nozzle - NPR = 5.0 - 32x27 Mesh - First-Order. 


00 + 30000 * 



NPR = 3.0 - 32x27 Mesh 



Figure 4.12a Velocity Vectors for Fully Deployed Thrust Reversing 
Nozzle - NPR = 3.0 - 32x27 Mesh - First-Order. 








Figure 4.14a Velocity Vectors for Fully Deployed Thrust Reversing 

Nozzle NPR = 3.0 - 32x27 Internal Mesh - 9x56 External 
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Figure 4.16 Variation of Normalized Thrust with Nozzle Pressure Ratio 
for a Fully Deployed Thrust Reversing Nozzle. 
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Figure 4.17 Sidewall Pressure Distribution at ZP = 0.203 cm - NPR 
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Figure 4.18 Sidewall Pressure Distribution at ZP = 0.838 cm - NPR 
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Figure 4.19 Sidewall Pressure Distribution at ZP - 1.473 cm - NPR 
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Figure 4.20 Pressures Along Flap for NPR 
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Figure 4.21 Mesh for 50% Deployed Thrust Reversing Nozzle 
22x20 Internal Mesh and 18x55 External Mesh 
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Figure 4.22 Pressure Contours for 50% Deployed Nozzle NPR= 
22x20 Internal Mesh and 18x55 External Mesh. 
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Velocity Vectors for Symmetric 50% Deployed Nozzle 
22x20 Internal Mesh and 18x55 External Mesh External 
Mach Number of Zero (Only External Vectors Shown) , 
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Mesh for 50% Deployed Thrust Reversing Nozzle with 
Downward Thrust Vectoring - 22x38 Internal Mesh an 
18x108 External Mesh. 
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50% Deployed Nozzle NPR= 
18x108 External Mesh. 
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l 0-30000 
00+30000 


148 



Reversing 




F igure 4 . 29 


Mesh for Steady State Flow Field Calculation 
Prior to Thrust Vectoring Angle Change. 


Time = 1.36 msec. 

Thrust Vector Angle = 0.0 deg. 
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F igur e 4 . 30 


Steady State Velocity Vectors Prior to Thrust 
Vectoring Angle Change. 


Time 

Thrust Vector Angle 


1.36 msec 
0.0 deg. 





Figure 4.32 Mesh for Transient Flow Field Calculation 
During Thrust Vectoring Angle Change. 


Time 

Thrust Vector Angle 


1.66 msec . 
28.90 deg. 




1 — s 


a — -a. 



Figure 4.33 Transient Flow Field Velocity Vectors During 
Thrust Vectoring Angle Change. 


Time 

Thrust Vector Angle 


1.66 msec . 
28.90 deg. 
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Figure 4.35 Mesh for Transient Flow Field Calculation 
After Thrust Vectoring Angle Change. 

Time = 1.71 msec. 

Thrust Vector Angle = 30.00 deg. 
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Figure 4.36 Transient Flow Field Velocity Vectors After 
I Thrust Vectoring Angle Change. 


Time 

Thrust Vector Angle 


1.71 msec . 
30.00 deg. 


/// 




Figure 4.37 Transient Flow Field P/P t Contours After 
Thrust Vectoring Angle Change. 

Time _ = 1.71 msec. 

Thrust Vectoring Angle = 30.00 deg. 




Figure 4.38 Mesh for Steady State Flow Field Calculation 
Prior to Transition from Partial to Full 
Thrust Reverser Deployment. 

Time = 1.44 msec . 

Rear Port Width = 0.0500 ft. 
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Figure 4.39 Steady State Flow Field Velocity Vectors Prior 
to Transition from Partial to Full Thrust 
Reverser Deployment. 


Time 

Rear Port Width 


1.44 msec 
0.0500 ft. 
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Figure 4.40 Steady State P/P*. Contours Prior to Transition 
from Partial to Full Thrust Reverser 
Deployment. 

Time = 1.44 msec. 

Rear Port Width = 0.0500 ft. 


i 




Figure 4.41 Mesh for Transient Flow Field Calculation 

^ During Transition from Partial to Full Thrust 

Reverser Deployment. 


Time 

Rear Port Width 


1.70 msec. 
0.0048 ft. 
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Figure 4.42a Transient Flow Field Velocity Vectors During 
Transition from Partial to Full Thrust 
Reverser Deployment 

Time = 1.70 msec. 

Rear Port Width = 0.0048 ft. 



Figure 4.42b Detail of Transient Flow Field Velocity 

Vectors Near Entrance to Rear Port (Boxed 
Region on Previous Plot). 
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Rear Port Width 


1.70 msec . 
0.0048 ft. 
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Figure 4.43 Transient Flow Field P/P t Contours During 
Transition from Partial to Full Thrust 
Reverser Deployment. 


Time 

Rear Port Width 


1.70 msec . 
0.0048 ft. 




Figure 4.44 Mesh for Transient Flow Field Calculation 
After Transition from Partial to Full 
Thrust Reverser Deployment. 


Time 

Rear Port Width 


1.79 msec 
0.0020 ft. 
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Figure 4.45 Transient Flow Field Velocity Vectors After 
Transition from Partial to Full Thrust 
Reverser Deployment. 

Time = 1.79 msec. 

Rear Port Width = 0.0020 ft. 
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Figure 4.46 Transient Flow Field P/P t Contours During 
Transition from Partial to Full Thrust 
Reverser Deployment. 
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APPENDIX A 

Diagonal izat ion of the Jacobians for the Inviscid Fluxes 
The inviscid fluxes through a surface of arbitrary 
orientation are 

/ " 

Pj ' S = I puq + P s x 

\ pvq + P s y 

( e+p) q (a.l) 

where q = us + vs and s and s are the components of a 
x y x y 

unit vector normal to the surface. 

The flux in equation a.l is a homogeneous function of 
the conservative variables U. So 
Pj - S = AU 

where A is the Jacobian of P^-S with respect to U. 
Following Warming et.al., (Reference 14) the Jacobian 
matrix, A, may be diagonalized by a set of three similarity 
transformat ions . 

A = T -1 R - 1 S _1 AS R T |s| 

Here T is a matrix which transforms incremental changes in 
conservative variables to incremental changes in 
nonconservative variables, R rotates from the global 
reference frame to a local reference frame with axis normal 
and tangent to the surface, and S transforms from 
incremental changes in the nonconservative variables to 
incremental changes in the characteristic variables. The 
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diagonal matrix, 

A= diag [u', u', u'+c, u'-c] 
contains the eigenvalues of the Jacobian matrix, A. 

The transformation matrices, and their inverses are 
given below. 


— 



■“ 

1 

0 

0 

0 

-u/p 

1/P 

0 

0 

-v/p 

0 

1/p 

0 

3(7-1 ) 

-u (7-1 ) 

-v (7-1 ) 

(7-1) 





1 

0 

0 

0 

u 

P 

0 

0 

V 

0 

P 

0 

a 

pu 

pv 

1 

(7-£_ 
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-s 
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s 
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0 
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0 
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S = 


1 

0 

0 

0 


0 

1 

0 

0 


0 

0 

pc 

-pc 


-1/c 

0 

1 

1 


.-1 


1 

0 

0 

0 


0 

1 

0 

0 


l/2c 

0 

l/2pc 

1/2 


l/2c 

0 

-l/2pc 

1/2 


In the above matrices 

a = 0.5 (u 2 + v 2 ) 
and c is the speed of sound. 
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APPENDIX B 

Block Pentadiagonal Solver 

The block pentadiagonal system of linear algebraic 
equations. Figure B.l, is solved using a block LU 
decomposition. The matrix is factored into two matrices, a 
block lower tridiagonal and a block upper triangular with 
identity matrices on the diagonal as shown in Figure B.2. 

The LU decomposition and inversion of the L matrix is 
done in one sweep in the direction of reducing the j index. 

1) A., * A.. , P-- = A.fc.., 0.. = A.} E.,, 


2) B jl-1 


' r jl 

= A-} Cjl , Sjl - 
, -1 

aT* e , 
31 :i 

y 3 l = 

A ,7R . , 

Dl Dl 


jl-1' 

A ji-i ' A n-i " 

B jl-l r jl 

r j 1-1 

■ A n-x “ii-i 

- B jl-1 0 

l E jl- 

!■ ^1-X - iji- 

l tR jl-l 


3) 

e . = b . - 
^3 3 

r j + 2 

4) 

A . = A . - 
3 3 

r j + l - 

5) 

r. = a' 1 

3 3 

lc : - b 3®s*i 1 ' «j ■ A 'j E r 


h - A ': lR j - Vj*i - D j y i+2 1 

6) repeat steps 3-5 for j = j 1 — 3 to 1 

In the above equation the vector [y] is equal to [L] ^ [ R ] . 
All that remains is to solve the equation 
[U] [ U] = [yl 

by sweeping in the direction of increasing j index. 

7) 5U 1 = y 1 
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8) 

5D 2 - *2 ' r 2 5U 1 


9) 

SU j ' - r 3 5U j-1 

- 5U j-2 

10) 

repeat step 9 for 

j = 3 to jl 


This algorithm requires the inversion of a 4x4 at each j 
location. This is done using a Crout decomposition 
(Reference 15 ) . 





Figure B.2 LU Decomposition of the Block Pentadiagonal Mat 
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APPENDIX C 


Ideal Mass Flow and Thrust 

The formula for the ideal mass flow is derived from 
★ ★ ★ 

m = uA = a A (c.l) 

where the * indicates the value of the quantity at a sonic 
throat. Using the isentropic relations 


Tjr- = 1 + M 2 and 



[1 



7 

7-1 


( c . 2 ) 

(c . 3 ) 


the definition of the speed of sound 

a 2 = 7RT , ( c . 4 ) 

the equation of state 

P = PRT , (c . 5 ) 

and the fact that the Mach number is unity at a sonic throat 
gives 


m . 

l 




1 

7-1 


r 2X_ RT l 1 / 2 

7 + l R T t 


* 

A 


( c . 6 ) 


The formula for the ideal thrust is obtained from 

F = ftue (c . 7 ) 

where ue is the flow at the exit of a ideally expanded 

nozzle (p = p ). The formula for u is obtained from 
r e oo e 


2 2 2 

u = a M 
e e e 


( c . 8 ) 


and equations c.2 - c.5. 


The result is 
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27, 


7-1 

e >7 


U e " { 7^t RT t 11 - <pT>' 


1/2 


(c . 9 ) 


Before substituting u g from equation c.9 into equation 
c.7 it should be decided which mass flow, A, to use. In the 
experimental results of Re and leavitt (Reference 1) and 
Carson, et.al., (Reference 2) the measured mass flow is used 
in the calculation of the ideal thrust. To be consistent, 
we will also use the measured mass flow when making 
comparisons with experiment. 

F = m T RT [1 - & 

1 Tit P t (C.10) 
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